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Whenever  band-limited  signals  are  measured  or  generated,  different 
distributions  of  a  fixed  number  of  receivers  and  transducers  lead  to  very 
different  resolutions.  A  Closely  related  set  of  issues  is  encountered  in 
the  numerical  solution  of  scattering  problems:  given  a  scatterer,  one 
would  like  to  find  nodes  on  its  surface  leading  to  roost  efficient  discret-j 
izations.  In  this  project,  we  have  constructed  numerical  algorithms  for 
the  determination  of  such  discretizations  in  one,  two,  and  three  dimen¬ 
sions,  and  used  the  obtained  results  for  the  design  of  optimal  phased 
antenna  arrays.  The  work  will  be  reported  in  four  papers;  copies  of  two 
of  these  papers  are  attached,  and  two  are  in  preparation.  Also,  a 
preliminary  patent  application  has  been  submitted. 
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1  Report 


Whenever  band-limited  signals  are  measured  or  generated,  the  locations  of  receivers 
or  transducers  have  to  be  selected;  it  is  well-known  that  different  distributions  lead  to 
very  different  resolutions  given  a  fixed  number  of  receivers  or  transducers.  A  closely 
related  set  of  issues  is  encountered  in  the  numerical  solution  of  scattering  problems: 
given  a  scatterer,  one  would  like  to  find  nodes  on  its  surface  leading  to  most  efficient 
discretizations. 

During  Phase  I  of  the  STTR  project  F49620-97-C-0052,  we  discovered  (somewhat 
serendipitously)  that  whenever  band-limited  signals  are  to  be  discretized,  measured,  or 
generated,  the  construction  of  optimal  (in  a  very  strong  sense)  configurations  of  nodes 
is  a  tractable  problem.  When  the  nodes  are  to  be  located  on  a  line  or  on  a  disk  in 
R^,  the  solution  is  a  fairly  straightforward  consequence  of  classical  results  obtained  by 
Slepian  and  his  collaborators  more  than  30  years  ago.  We  have  constructed  the  necessary 
numerical  tools,  which  are  quite  efficient;  the  resulting  discretizations  are  a  dramatic 
improvement  over  the  ones  currently  employed. 

The  basic  analytical  apparatus  for  dealing  with  band-limited  functions  are  the  clas¬ 
sical  Prolate  Spheroidal  Wave  Functions  (PSWFs)  and  their  generalizations.  We  were 
surprized  to  discover  that  the  existing  numerical  tools  for  the  evaluation  of  PSWFs  leave 
much  to  be  desired,  being  based  on  the  so-called  Bouwkamp  algorithm,  constructed 
in  1941  (see  [1]).  Indeed,  a  fairly  straightforward  analysis  shows  that  the  Bouwkamp 
scheme  is  numerically  unstable,  except  for  small-scale  problems;  as  a  result,  there  ap¬ 
pears  to  exist  a  belief  (see,  for  example,  [18])  that  the  numerical  evaluation  of  PSWFs 
presents  severe  numerical  difficulties.  When  we  examined  the  Bouwkamp  algorithm,  we 
discovered  that  a  very  simple  alteration  eliminates  the  instability  completely.  In  fact, 
the  required  change  is  no  more  than  the  use  of  standard  modern  numerical  techniques  for 
the  solution  of  a  fairly  simple  Sturm-Liouville  problem.  Needless  to  say,  such  techniques 
did  not  exist  in  1941,  when  C.  J.  Bouwkamp  developed  his  scheme.  In  the  process  of  de¬ 
signing  the  requisite  numerical  tools,  we  discovered  that  the  PSWFs  posess  an  extremely 


rich  collection  of  analytical  properties;  we  list  some  of  these  properties  in  [20],  where  we 
also  describe  interpolation  and  quadrature  techniques  for  band-limited  functions. 

Next,  we  applied  the  constructed  apparatus  to  the  design  of  configurations  of  trans¬ 
ducers  for  linear  phased  antenna  arrays,  and  to  antenna  arrays  on  disks  in  the  plane.  This 
work  is  reported  in  [19]  (attached).  Construction  of  optimal  configurations  of  nodes  on 
more  complicated  regions  requires  additional  mathematical  apparatus;  such  apparatus 
has  been  largely  designed,  and  largely  implemented  numerically.  The  paper  describing 
it  is  in  preparation.  We  are  also  investigating  applications  of  the  constructed  numerical 
techniques  in  the  design  of  receiver  configurations  in  seismic  data  collection  (as  encoun¬ 
tered  in  oil  exploration  and  related  areas),  electronic  beam  steering,  and  several  other 
environments.  A  preliminary  patent  application  has  been  filed  by  the  Office  of  Cooper¬ 
ative  Research  at  Yale  University. 

2  Appendix:  Antenna  Patterns  and  Corresponding 
Optimal  Element  Distributions 

2.1  Characteristics  of  an  antenna  pattern 

Depending  on  the  situation,  the  design  of  an  antenna  array  attempts  to  optimize  certain 
characteristics  of  the  resulting  far-field  pattern,  subject  to  certain  constraints  on  the 
number,  power,  etc.  of  the  elements.  Since  the  principal  purpose  of  this  work  is  to 
develop  a  technique  for  the  selection  of  the  locations  of  the  elements  that  approximate  a 
user-specified  pattern,  we  could  use  any  reasonable  far-field  pattern  to  be  approximated. 
In  subsection  2.2,  2.3,  we  construct  optimal  element  distributions  for  the  so-called  sector 
patterns  and  cosecant  pattern,  respectively;  a  detailed  discussion  of  these  (and  several 
other)  pattern  cans  be  found,  for  example  in  [6]. 

We  will  say  that  the  antenna  pattern  has  the  e-bandwidth  b  if 
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in  other  words,  the  proportion  of  the  energy  radiated  outside  the  e-beamwidth  from  the 
axis  of  the  beam  is  equal  to  e.  The  supergain  of  an  antenna  is  defined  (see,  for  example, 
[22]),  as  the  ratio 
+  00 

/  \F'{x)\'^  dx 

- .  (2) 

f  \F{x)\^  dx 
-1 

The  supergain  (sometimes  referred  to  as  superdirectivity)  measures  the  ratio  of  the  en¬ 
ergy  associated  with  the  total  spectrum  of  the  antenna  to  the  energy  in  its  visible  spec¬ 
trum;  while  detailed  discussion  of  supergain  and  related  issues  is  outside  the  scope  of  this 
report,  we  will  observe  that  antenna  arrays  with  large  degrees  of  supergain  would  violate 
the  uncertainty  principle,  and  thus  are  physically  impossible.  Attempts  to  construct 
supergain  antennae  result  in  rapidly  (exponentially)  growing  Ohmic  losses,  prohibitive 
accuracy  requirements,  extremely  low  bandwidth,  etc.  Thus,  any  potentially  useful  pro¬ 
cedure  for  the  design  of  antenna  arrays  has  to  limit  the  supergain  of  the  resulting  patterns. 

2.2  Sector  patterns 

It  is  often  desirable  to  construct  antenna  patterns  that  are  as  constant  as  possible  within 
the  main  beam,  and  as  small  as  possible  outside  it;  in  other  words,  ideally,  the  pattern 


would  be  defined  by  the  formulae 

Fb{x)  =  1  for  \x\  <  b, 

(3) 

Fb{x)  =  0  for  |a;|  >  6, 

(4) 

with  b  a  real  number  such  that  0  <  b  <  k.  Needless  to  say,  the  function  Ft  defined  by 
the  formulae  (3),  (4)  is  not  band-limited,  and  some  approximation  has  to  be  used.  A 
standard  procedure  is  to  truncate  the  Fourier  Transform  of  Ft,  approximating  it  by  the 
function  F  defined  by  the  formula 
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(see,  for  example,  [21]).  An  important  special  case  occurs  when  b  =  k,  with  (5)  assuming 
the  form 


Fk(x 


smjk  ■  t)  _ 


(6) 


obviously,  the  latter  expression  is  a  band-limited  approximation  of  the  5-function.  An¬ 
other  frequently  encountered  situation  is  that  of  6  =  k/2,  so  that  (5)  assumes  the  form 


sm(|  •  t) 


Akx-t 


(7) 


which  is  a  band-limited  approximation  to  the  beam  that  is  equal  to  1  for  — 1/2  <  a:  <  1/2 
and  to  zero  elsewhere. 

In  Section  2.4  below,  we  demonstrate  optimal  element  configurations  that  produce 
approximations  to  the  patterns  (6),  (7)  with  k  =  207r,107r,  32.46  767r. 


Remark  2.1  While  (5)  is  by  no  means  the  only  possible  band-limited  approximations  to 
to  Ffc,  it  is  quite  satisfactory  in  most  cases,  in  addition  to  being  simple.  Furthermore,  the 
principal  purpose  of  this  report  is  to  describe  a  technique  for  the  selection  of  locations  of 
the  nodes,  given  a  pattern  to  be  approximated.  Thus,  we  ignore  the  issue  of  the  optimal 
choice  of  Fi,. 


2.3  Cosecant  patterns 

Another  standard  far-field  radiation  pattern  is  the  so-called  cosecant  pattern  (see,  for 
example,  [8]).  Given  two  real  numbers  0  <  a  <  6  <  1,  the  cosecant  pattern  Fa,b  is  defined 
by  the  formula 

FaM  =  I  (8) 

for  all  X  G  [a,  b],  and 

Fa,b{x)  =  0  (9) 
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for  all  X  e  ([-1,1]  \  [a,  6]).  Again,  the  function  Fa^b  defined  by  the  formulae  (8),  (9) 
is  not  band-limited.  Before  our  scheme  can  be  applied  to  the  latter  has  to  be 
approximated  with  a  band-limited  function;  as  discussed  in  Section  2.1  above,  if  such 
an  approximation  is  to  be  useful  as  an  antenna  pattern,  its  supergain  factor  has  to  be 
controlled.  Fortunately,  a  procedure  for  such  an  approximation  has  been  in  existence  for 
more  than  35  years  (see,  [7]);  the  algorithm  of  [7]  is  a  modification  of  the  least-squares  ap¬ 
proach  permitting  the  user  to  limit  the  supergain  factor  of  the  obtained  pattern  explicitly. 
At  the  time,  the  utility  of  the  scheme  of  [7]  was  limited  by  the  (perceived)  difficulty  in 
the  numerical  evaluation  of  Prolate  Spheroidal  Wave  functions;  given  the  present  state 
of  numerical  analysis,  this  difficulty  is  non-existent,  and  it  is  this  author’s  impression 
that  the  insights  of  [7],  [8]  deserve  more  attention  than  they  have  been  receiving. 

2.4  Optimal  distributions  of  elements 

In  this  subsection,  we  briefly  describe  an  algorithm  for  the  construction  of  optimal  (in 
the  sense  defined  below)  element  configurations  for  the  generation  of  antenna  patterns, 
of  which  the  patterns  (4)- (6)  are  special  cases.  As  will  be  seen,  the  procedure  is  in  fact 
applicable  to  the  design  of  element  configurations  for  very  general  far-field  patterns. 

We  start  with  observing  that  the  far-field  pattern  F  is  an  integral  over  the  interval 
[—1,1]  of  functions  of  the  form 

(J{u)  ■  (10) 

with  X  =  cos{0)  determined  by  the  direction  9  in  which  the  far-field  is  being  evaluated.  In 
other  words,  the  problem  of  finding  efficient  antenna  element  distributions  is  equivalent  to 
that  of  constructing  quadrature  formulae  for  functions  of  the  form  (10).  In  our  numerical 
experiments,  the  techniques  of  [2])  (after  some  tuning)  have  always  been  successful  in 
finding  the  Gaussian  quadratures  for  integrals  of  the  form  (3);  some  of  our  results  are 
presented  in  Section  3  below. 


3  Numerical  Examples 


In  this  section,  we  present  examples  of  optimal  element  distributions  generating  the 
patterns  of  the  preceding  Section;  all  of  the  results  presented  here  have  been  obtained 
numerically.  Antenna  patterns  we  present  are  compared  to  the  antenna  patterns  given 
by  uniform  source  distributions;  configurations  of  elements  approximating  these  antenna 
patterns  are  compared  to  equispaced  distributions  of  elements  generating  the  same  an¬ 
tenna  patterns. 

3.1  Optimal  distributions  of  elements 

In  this  section,  we  demonstrate  the  results  of  the  application  of  the  techniques  of  Sec¬ 
tion  2.4  of  this  report  to  the  types  of  antenna  patterns  described  in  the  Sections  2.2,  2.3. 

In  all  cases,  we  choose  the  size  of  an  antenna  array  and  a  pattern  to  be  reproduced,  and 
use  the  scheme  outlined  in  Section  2.4  to  design  a  distribution  of  antenna  elements  (both 
the  locations  and  the  intensities)  located  within  the  chosen  array  that  reproduces  the 
required  pattern.  For  comparison,  we  also  generate  optimal  (in  the  least  squares  sense) 
approximations  to  the  desired  pattern  generated  by  equispaced  elements  located  within 
the  same  array.  Since  the  number  of  equispaced  nodes  required  to  obtain  a  reasonable 
approximation  to  the  desired  pattern  is  (in  many  cases)  much  greater  than  the  number  of 
optimally  chosen  nodes,  for  each  example  we  demonstrate  patterns  generated  by  several 
such  configurations.  In  this  manner,  the  numbers  of  optimally  chosen  nodes  necessary 
to  obtain  reasonable  approximations  to  the  desired  patterns  can  be  compared  to  the 
numbers  of  equispaced  nodes  required  to  obtain  similar  results. 

3.1.1  Sector  patterns 

Example  3.1  The  first  example  we  consider  is  of  the  pattern  defined  by  the  formula  (7), 
with  k  =  62.8312,  so  that  the  size  of  the  array  is  20  wavelengths. 

In  Figure  5,  we  display  an  approximation  to  the  pattern  obtained  with  19  elements, 
overlayed  with  the  exact  pattern;  the  locations  of  the  elements  are  displayed  in  Figure  5a; 
the  relative  error  of  the  obtained  approximation  is  5.01%. 


Similarly,  in  Figure  5g,  we  display  the  approximation  to  the  pattern  obtained  with  21 
elements,  overlayed  with  the  exact  pattern;  the  relative  error  of  the  obtained  approxima¬ 
tion  is  0.443%;  in  Figure  5h,  we  display  the  the  approximation  obtained  with  17  elements. 
In  the  latter  case,  the  relative  error  of  the  obtained  approximation  is  6.43%;  Figure  5i 
depicts  the  17 -node  distribution  producing  the  approximation  illustrated  in  Figure  5h. 
Finally,  Figure  5j  contains  a  graph  of  the  values  of  the  sources  located  at  the  17  nodes 
depicted  in  Figure  5i  and  generating  the  pattern  shown  in  Figure  5h. 

For  comparison,  the  optimal  approximation  obtained  with  19,  24,  29,  31,  and  34 
equispaced  elements  are  displayed  in  Figures  5b,  5c,  5d,  5e,  5f,  respectively ;  these  are 
also  overlayed  with  the  exact  pattern. 

Example  3.2  Our  second  example  is  identical  to  the  first  one,  with  the  exception  that 
k  =  31.416,  so  that  the  size  of  the  array  is  10  wavelengths. 

In  Figure  6,  we  display  an  approximation  to  the  pattern  obtained  with  9  elements, 
overlayed  with  the  exact  pattern;  the  locations  of  the  elements  are  displayed  in  Figure  6a; 
the  relative  error  of  the  obtained  approximation  is  11.2%. 

Similarly,  in  Figure  6f,  we  display  the  approximation  to  the  pattern  obtained  with  11 
elements,  overlayed  with  the  exact  pattern;  the  relative  error  of  the  obtained  approxima¬ 
tion  is  0.600%. 

For  comparison,  the  optimal  approximation  obtained  with  9,  14,  16,  and  18  pquispaced 
elements  are  displayed  in  Figures  6b,  6c,  6d,  5e,  respectively;  these  are  also  overlayed 
with  the  exact  pattern. 

Example  3.3  Our  third  example  is  identical  to  the  preceding  two,  with  the  exception 
that  k  =  102,  so  that  the  size  of  the  array  is  about  32.45  wavelengths. 

In  Figure  7a,  we  display  an  approximation  to  the  pattern  obtained  with  23  optimally 
distributed  elements,  overlayed  with  the  exact  pattern  and  with  the  pattern  obtained  with 
23  equispaced  elements. 

The  relative  error  of  the  obtained  approximation  is  5.4%;  needless  to  say,  the  error  of 
the  approximation  obtained  with  the  equispaced  nodes  is  more  than  70%.  ^45  can  be  seen 


from  Figure  7c,  the  actual  size  of  the  obtained  23-element  array  is  about  21  wavelengths; 
in  other  words,  in  order  to  obtain  this  precision,  the  array  needs  to  be  about  2/3  of  the 
nominal  (maximum  permitted)  length. 

In  Figure  7b,  we  display  the  approximation  to  the  pattern  obtained  with  42  and  48 
elements,  overlayed  with  the  exact  pattern. 

It  is  worth  noting  that  with  33  optimally  distributed  elements,  the  pattern  is  approxi¬ 
mated  to  the  precision  0.12%;  we  do  not  display  the  obtained  pattern  since  it  is  visually 
indistinguishable  from  the  pattern  being  approximated. 

Example  3.4  Our  final  example  is  somewhat  different  from  the  preceding  ones,  in  that 
instead  of  approximating  a  sector  pattern,  we  approximate  a  cosecant  pattern  (see  (8),  (9) 
in  Subsection  2.3  above). 


In  this  example,  we  set 

a  =  sm(15°), 

(11) 

b  =  sin(75°), 

(12) 

and  use  the  procedure  of  [7]  to  approximate  Fa,b  with  a  band-limited  function.  The  band- 
limit  has  been  more  or  less  arbitrarily  set  to  110,  resulting  in  an  antenna  array  about  35 
wavelengths  in  size,  and  the  supergain  factor  of  the  approximation  was  set  to  1.1. 

In  Figure  8a,  we  display  an  approximation  to  the  pattern  obtained  with  53  optimally 
distributed  elements,  overlayed  with  the  exact  bandlimited  pattern  and  with  the  pattern 
obtained  with  53  equispaced  elements. 

The  relative  error  of  the  obtained  approximation  is  1.79%;  the  error  of  the  approxi¬ 
mation  obtained  with  the  equispaced  nodes  is  about  42%. 

In  Figure  8b,  we  display  the  approximation  to  the  pattern  obtained  with  47  optimally 
distributed  elements,  overlayed  with  the  exact  pattern;  the  purpose  of  this  final  figure  is 
to  demonstrate  the  behavior  of  the  scheme  when  the  number  of  elements  is  insufficient 
(i.e.  when  the  array  is  underresolved). 

It  is  worth  noting  that  it  takes  about  70  equispaced  nodes  to  obtain  the  resolution 
obtained  with  47  optimally  chosen  ones. 


The  following  observations  can  be  made  from  Figures  5  -  8b,  and  from  the  more 
detailed  numerical  experiments  performed  by  the  author. 

1.  In  order  to  obtain  reasonable  precision,  the  scheme  requires  about  1  point  per  wave¬ 
length  in  the  antenna  array;  this  is  more  or  less  independent  from  the  structure  of  the 
beam  as  long  as  the  pattern  is  symmetric  about  the  point  x  =  0.  This  fact  is  observed 
numerically,  even  for  modest  numbers  of  nodes;  for  large-scale  arrays,  this  statement 
(interpreted  asymptotically)  can  be  proved  rigorously.  For  certain  beam  structures,  the 
required  number  of  nodes  is  even  less  (see  Example  3.3).  The  reasons  for  these  additional 
savings  are  subtle,  and  have  to  do  with  the  fact  that  the  continuous  source  distribution 
generating  the  pattern  is  relatively  small  on  a  large  part  of  the  antenna  array;  the  al¬ 
gorithm  of  [2]  takes  advantage  of  this  fact  to  reduce  the  number  of  nodes.  When  the 
beam  is  not  symmetric  about  x  =  0,  the  number  of  elements  required  does  depend  on 
the  structure  of  the  pattern,  and  the  dependence  is  fairly  complicated.  Generally,  the 
improvement  for  non-symmetric  beams  is  less  than  that  for  the  symmetric  ones. 

2.  The  qualiative  behavior  of  the  scheme  is  similar  to  that  of  the  Gaussian  quadratures 
in  that  it  displays  no  convergence  at  all  until  a  certain  minimum  number  of  nodes  is 
achieved;  after  that,  the  convergence  is  very  fast.  This  behavior  is  not  surprising,  since 
the  scheme  is  based  on  a  Generalized  Gaussian  quadrature. 

3.  For  the  sector  pattern  with  the  sector  [—1/2, 1/2],  the  scheme  reduces  the  required 
number  of  nodes  by  a  factor  of  about  1.5  for  small-scale  problems,  and  roughly  by  a 
factor  of  2  for  large-scale  ones;  again,  for  large-scale  problems,  an  asymptotic  version  of 
this  statement  can  be  proven  rigorously. 

4.  For  the  cosecant  pattern  with  the  parameters  specified  by  (11),  (12),  the  number 
of  nodes  required  is  reduced  by  approximately  a  factor  of  1.4.  As  the  sidelobe  level  is 
reduced,  the  improvement  obtained  by  going  from  the  equispaced  discretization  to  the 
optimal  one  increases  rapidly. 

5.  An  examination  of  Figures  5a,  6a  shows  that  while  the  optimal  nodes  are  by  no  means 
uniform,  they  display  no  clustering  behavior. 


6.  An  examination  of  Figure  5j  shows  that  the  intensities  of  individual  elements  do  not 
become  large;  this  is  confirmed  by  the  more  extensive  numerical  experiments  performed 
by  the  author. 

7.  The  combination  of  the  preceding  two  paragraphs  (combined  with  additional  numer¬ 
ical  experiments  and  analysis)  provide  evidence  that  configurations  of  this  type  should 
pose  no  supergain  problems. 

4  Generalizations 

The  results  described  above  admit  radical  generalizations  in  several  directions;  several 
such  directions  are  discussed  below, 

1.  Conformal  one-dimensional  arrays.  The  extension  of  the  techniques  of  this  report 
to  one-dimensional  arrays  located  on  curves  in  is  completely  straightforward,  involving 
only  a  modest  increase  of  the  CPU  time  requirements  of  the  procedure.  Improvement  in 
the  number  of  nodes  required  to  produce  a  prescribed  pattern  is  similar  to  that  in  the 
case  of  a  linear  array. 

2.  Planar  two-dimensional  arrays.  A  straightforward  generalization  of  the  results  of 
Sections  2,  3,  is  to  rectangular  planar  arrays.  Here,  a  tensor  product  quadrature  can  be 
constructed  from  the  quadratures  of  Sections  2,  3,  possessing  all  of  the  desirable  prop¬ 
erties  of  the  latter.  Obviously,  the  advantage  in  the  number  of  transducers  is  squared, 
so  that  (for  example)  replacing  50  nodes  in  each  of  the  two  directions  by  23  nodes  (see 
Example  3.3  above)  will  lead  to  a  factor  of  (50/23)^  ~  4.7  savings  in  the  number  of 
elements. 

The  theory  of  Section  2  has  been  extended  for  disk-shaped  arrays,  via  {inter  alia)  the 
techniques  developed  in  [12].  The  improvement  in  the  number  of  nodes  is  comparable  to 
that  obtained  in  the  rectangular  geometry,  and  the  CPU  time  requirements  do  not  differ 
appreciably  from  those  in  the  case  of  linear  one-dimensional  arrays. 


The  extension  of  the  theory  to  more  general  geometries  in  the  plane  is  in  progress.  At 
the  present  time,  our  only  numerical  experiments  have  been  with  arrays  on  triangles;  the 
results  are  encouraging,  but  the  CPU  time  requirements  of  the  algorithms  are  excessive 
(we  have  only  been  able  to  design  triangular  arrays  about  6  wavelengths  in  size).  We 
are  now  in  the  process  of  constructing  a  more  efficient  numerical  procedure  for  such 
computations. 

3.  Conformal  two-dimensional  arrays.  The  only  environment  in  which  we  have 
a  satisfactory  theory  is  when  the  array  is  located  on  a  surface  of  revolution;  even  in 
this  environment,  no  experiments  have  been  performed.  We  have  not  investigated  more 
general  conformal  two-dimensional  arrays  in  sufficient  detail. 


Figure  5;  The  pattern  created  by  the  19  optimal  elements,  depicted  in  Figure 
5a  as  described  in  Example  3.1 
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Figure  5a:  The  distribution  of  elements  creating  the  pattern  depicted  in 
Figure  5,  as  described  in  Example  3.1 


Figure  5b:  The  optimal  approximation  to  the  sector  pattern  generated  by  19 
equispaced  nodes,  as  described  in  Example  3.1 


Figure  5c:  The  optimal  approximation  to  the  sector  pattern  generated  by  24 
equispaced  nodes,  as  described  in  Example  3.1 


Figure  5d:  The  optimal  approximation  to  the  sector  pattern  generated  by  29 
equispaced  nodes,  as  described  in  Example  3.1 


Figure  5e:  The  optimal  approximation  to  the  sector  pattern  generated  by  31 
equispaced  nodes,  as  described  in  Example  3.1 


Figure  5f;  The  optimal  approximation  to  the  sector  pattern  generated  by  34 
equispaced  nodes,  as  described  in  Example  3.1 


Figure  5g:  The  optimal  approximation  to  the  sector  pattern  generated  by  21 
optimal  nodes,  as  described  in  Example  3.1 


Figure  5h:  The  optimal  approximation  to  the  sector  pattern  generated  by  17 
optimal  nodes,  as  described  in  Example  3.1 
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Figure  5i:  The  distribution  of  17  elements  creating  the  pattern  depicted  in 
Figure  5h,  as  described  in  Example  3.1 


1.4 


Figure  5j :  The  values  of  the  sources  located  at  the  nodes  depicted  in  Figure  5i 
and  generating  the  pattern  depicted  in  Figure  5h,  as  described  in  Example  3.1 
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Figure  6b;  The  optimal  approximation  to  the  sector  pattern  generated  by  9 
equispaced  nodes,  as  described  in  Example  3.2 


Figure  6c:  The  optimal  approximation  to  the  sector  pattern  generated  by  14 
equispaced  nodes,  as  described  in  Example  3.2 


Figure  6d;  The  optimal  approximation  to  the  sector  pattern  generated  by  16 
equispaced  nodes,  as  described  in  Example  3.2 


Figure  6e:  The  optimal  approximation  to  the  sector  pattern  generated  by  18 
equispaced  nodes,  as  described  in  Example  3.2 
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Figure  7a:  The  approximation  to  the  sector  pattern  generated  by  23  optimal 
elements,  vs.  optimal  approximation  by  23  equispaced  nodes,  as  described  in 

Example  3.3 


Figure  7b:  The  optimal  approximations  to  the  sector  pattern  generated  by  42 
and  48  equispaced  nodes,  as  described  in  Example  3.3 
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Whenever  physical  signals  are  measured  or  generated,  the  locations  of  receivers  or  trans¬ 
ducers  have  to  be  selected.  Most  of  the  time,  this  appears  to  be  done  on  an  ad  hoc  basis. 
For  example,  when  a  string  of  geophones  is  used  in  the  measurements  of  seismic  data  in 
oil  exploration,  the  receivers  are  located  at  equispaced  points  on  an  interval.  When  phased 
array  antennae  are  constructed,  their  shapes  are  determined  by  certain  aperture  consid¬ 
erations;  round  and  rectangular  shapes  are  common.  When  antenna  beams  are  steered 
electronically,  it  is  done  by  changing  the  phases  (and  sometimes,  the  amplitudes)  of  the 
transducers.  Again,  these  transducers  are  located  in  a  region  of  predetermined  geometry, 
and  their  actual  locations  within  that  geometry  are  chosen  via  some  heuristic  procedure. 
In  all  these  (and  many  other)  cases,  the  signals  being  received  or  generated  are  band-limited. 
Optimal  representation  of  such  signals  has  been  studied  in  detail  by  Slepian  et.  al.  more 
than  30  years  ago,  and  some  of  the  obtained  results  were  applied  by  D.  Rhodes  to  the 
design  of  antenna  patterns;  further  development  of  this  line  of  research  appears  to  have 
been  hindered  by  the  absence  at  the  time  of  necessary  numerical  tools.  We  combine  these 
classical  results  with  the  recently  developed  apparatus  of  Generalized  Gaussian  Quadratures 
to  construct  optimal  nodes  for  the  measurement  and  generation  of  band-limited  signals.  In 
this  report,  we  describe  the  procedure  based  on  these  techniques  for  the  design  of  such 
receiver  (and  transducer)  configurations  in  a  variety  of  environments. 
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1  Introduction 


When  measurements  are  performed,  it  often  happens  that  the  signal  to  be  measured  is 
well  approximated  by  linear  combinations  of  oscillatory  exponentials,  i.e.  functions  of 
the  form 

(1) 

j=i 

in  one  dimension,  of  the  form 

(2) 

j=i 

in  two  dimensions,  and  of  the  form 

(3) 

j=i 

in  three  dimensions.  In  most  cases,  the  signal  is  band-limited,  i.e.  there  exist  such  real 
positive  a  that  all  I  <  j  <  n, 

I  A,  1<  a  (4) 

in  one  dimension, 

Aj  +  (5) 

in  two  dimensions,  and 

+  +  (6) 

in  three  dimensions. 

As  is  well-known,  most  measurements  of  electromagnetic  and  acoustic  data  (espe¬ 
cially  at  reasonably  high  frequencies)  are  of  this  form.  Examples  of  such  situations 
include  geophone  and  hydrophone  strings  in  geophysics,  phased  array  antennae  in  radar 


systems,  multiple  transceivers  in  ultrasound  imaging,  and  a  number  of  other  applications 
in  astrophysics,  medical  imaging,  non-destructive  testing,  etc. 

In  this  report,  we  describe  a  procedure  for  determining  the  optimal  distribution  of 
sources  and  receivers  that  maximizes  accuracy  and  resolution  in  measuring  band-limited 
data  given  a  fixed  number  of  receivers.  Alternatively,  the  procedure  can  be  used  to 
determine  the  optimal  distribution  of  receivers  that  will  minimize  their  number  given 
specified  accuracy  and  resolution.  While  the  techniques  described  in  this  note  are  fairly 
general,  we  describe  them  in  detail  in  the  case  of  linear  antenna  arrays;  the  changes 
needed  to  generalize  the  approach  to  other  cases  are  summarized  in  Section  6. 

Remark  1.1  One  of  principal  issues  in  the  design  of  antenna  arrays  is  the  treatment 
(or  avoidance)  of  the  so-called  supergain  (or  superdirectivity).  Supergain  is  the  con¬ 
dition  that  occurs  when  an  antenna  design  is  attempted  that  is  prohibited  (or  nearly 
prohibited)  by  the  Heisenberg  principle;  technically,  it  occurs  in  the  form  of  very  closely 
spaced  elements  operating  out  of  phaze,  and  leads  to  prohibitive  Ohmic  losses  in  trans¬ 
mitting  antennae,  loss  of  sensitivity  in  receiving  ones,  etc.  Since  the  purpose  of  this 
note  is  to  introduce  techniques  for  selecting  the  locations  of  elements  for  a  prescribed 
antenna  pattern,  we  avoid  the  issue  of  choosing  the  antenna  pattern  altogether.  Instead, 
we  observe  design  optimal  element  distributions  for  several  standard  far-field  patterns 
(see  Section  5.1),  and  we  observe  that  the  scheme  for  choosing  optimal  distributions  of 
elements  is  virtually  independent  of  the  patterns  being  approximated. 

Technically,  the  approach  taken  here  is  to  observe  that  designing  an  antenna  array 
can  be  viewed  as  constructing  a  quadrature  formula  for  the  integration  of  certain  special 
classes  of  functions.  Using  recently  developed  techniques  for  the  construction  of  so-called 
Generalized  Gaussian  Quadratures,  we  obtain  both  nodes  and  weights  that  are  optimal 
(in  a  very  strong  sense)  for  the  required  antenna  pattern. 

The  structure  of  this  note  is  as  follows.  In  Section  2,  we  summarize  some  of  the  math¬ 
ematical  apparatus  to  be  used:  Ghebychev  Systems,  Generalized  Gaussian  Quadratures, 
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etc.  In  Section  3,  we  recapitulate  some  of  the  standard  antenna  theory,  primarily  to 
introduce  the  necessary  notation.  In  Section  4,  element  distributions  given  a  specific  an¬ 
tenna  pattern.  In  Section  5,  we  illustrate  our  approach  with  several  numerical  examples, 
and  Section  6  contains  a  discussion  of  the  generality  of  the  schemes  presented. 

2  Analytical  Preliminaries 

In  this  section,  we  summarize  several  known  facts  about  classical  Special  functions.  All 
of  these  facts  can  be  found  in  the  literature;  detailed  references  are  given  in  the  text. 


2.1  Chebyshev  systems 


Definition  2.1  A  sequence  of  functions  will  be  referred  to  as  a  Chebyshev 

sj'stem  on  the  interval  [a,  6]  if  each  of  them  is  continuous  and  the  determinant 


•••  M^n) 

•^n(^l)  '  ■  ■ 

is  nonzero  for  any  sequence  of  points  Xi,. ..  ,Xn  such  that  a  <  Xi  <  X2 . . .  <  Xn  <  b. 


(7) 


An  alternate  definition  of  a  Chebyshev  system  is  that  any  linear  combination  of  the 
functions  with  nonzero  coefficients  must  have  no  more  than  n  zeros. 

Examples  of  Chebyshev  and  extended  Chebyshev  systems  include  the  following  (ad¬ 
ditional  examples  can  be  found  in  [8]). 


Example  2.1  The  powers  l,x,x^, . . .  ,x'^  form  an  extended  Chebyshev  system  on  the 
interval  {—oo,  oo). 


Example  2.2  The  exponentials  e  e  . . . ,  e  form  an  extended  Chebyshev  sys¬ 
tem  for  any  Ai, . . . ,  >  0  on  the  interval  [0,  cxd). 


Example  2.3  The  functions  l,cosa;,  sina;,cos2a;,  sin2x, . . . ,  cos  nx,  sin  nx  form  a  Cheby¬ 
shev  system  on  the  interval  [0,27r]. 


Example  2.4  Suppose  that  c  >  0  is  a  real  number,  w  is  a  positive  function  [— 1, 1]  — >  II 
such  that  w  G  c^[— 1, 1]  and  w{—x)  —  w{x)  for  all  x  G  [—1, 1],  n  is  a  natural  number, 
and  the  operators  P,  Q  :  1, 1]  — >  1]  are  defined  by  the  formulae 

P{(l)){x)  ~  J  w{t)  •  (j){t)  dt  (8) 

Q  =  P*oP.  (9) 

Suppose  further  that  (f)i,(f)2,  ■  ■  ■  are  the  eigenfunctions  of  Q,  Ai,A2,...  are  the  corre¬ 
sponding  eigenvalues,  and  Ai  >  A2  >  A3....  Then  all  eigenfunctions  of  Q  (also  known 
as  the  right  singular  vectors  of  P)  can  be  chosen  to  be  real.  Furthermore,  the  functions 
■  ■,4>n  constitute  a  Chebychev  system  on  the  interval  [—1,1]. 

2.2  Generalized  Gaussian  quadratures 

A  quadrature  rule  is  an  expression  of  the  form 


(10) 


j=i 


where  the  points  xj  G  R.  and  coefficients  wj  G  R.  are  referred  to  as  the  nodes  and  weights 
of  the  quadrature,  respectively.  They  serve  as  approximations  to  integrals  of  the  form 

rb 

/  4>{x)  ■  ui{x)dx  (11) 

Ja 

with  u)  is  an  integrable  non-negative  function. 

Quadratures  are  typically  chosen  so  that  the  quadrature  (10)  is  equal  to  the  desired 
integral  (11)  for  some  set  of  functions,  commonly  polynomials  of  some  fixed  order.  Of 
these,  the  classical  Gaussian  quadrature  rules  consist  of  n  nodes  and  integrate  polynomi¬ 
als  of  order  2n  —  1  exactly.  In  [13],  the  notion  of  a  Gaussian  quadrature  was  generalized 
as  follows: 

Definition  2.2  A  quadrature  formula  will  be  referred  to  as  Gaussian  with  respect  to  a 
set  of  2n  functions  (f)i, ,  4>2n  '■  [a,  6]  — >■  R  and  a  weight  function  u>  ;  [a,  b]  R+,  if  it 
consists  of  n  weights  and  nodes,  and  integrates  the  functions  (fi  exactly  with  the  weight 
function  u>  for  alii  =  1, . . . ,  2n.  The  weights  and  nodes  of  a  Gaussian  quadrature  will  be 
referred  to  as  Gaussian  weights  and  nodes  respectively. 
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<c 


The  following  theorem  appears  to  be  due  to  Markov  [15,  16];  proofs  of  it  can  also  be 
found  in  [10]  and  [8]  (in  a  somewhat  different  form). 

Theorem  2.1  Suppose  that  the  functions  <^i,...,^2n  ^  R,  form  a  Chehyshev 

system  on  [a,  5] .  Suppose  in  addition  that  w  :  [a,  5]  — >■  R  is  a  non-negative  integrable 
function  [a,  b]  R.  Then  there  exists  a  unique  Gaussian  quadrature  for  the  functions 
(f>i,  ■  ■  ■  1  ^2n  on  [a,  b]  with  respect  to  the  weight  function  u.  The  weights  of  this  quadrature 
are  positive. 

Remark  2.1  While  the  existence  of  Generalized  Gaussian  Quadratures  was  observed 
more  than  100  years  ago,  the  constructions  found  in  [15,  16],  [3,  10],  [7,  8]  do  not  easily 
yield  numerical  algorithms  for  the  design  of  such  quadrature  formulae;  such  algorithms 
have  been  constructed  recently  (see  [13,  28,  2]).  The  version  of  the  procedure  found  in 
[2]  was  used  to  produce  the  results  presented  in  the  Examples  5.1,  5.2,  5.3  in  Section  5.1; 
the  reader  is  referred  to  [2]  for  details. 

Applying  Theorem  2.1  to  the  Example  2.4,  we  obtain  the  following  theorem. 

Theorem  2,2  Suppose  that  under  the  conditions  of  Example  2.4,  n  is  even.  Then 
there  exist  n/2  points  ti,t2, .  ■ .  ,tn/2  on  the  interval  [—1,1]  and  positive  real  numbers 
Wi,W2,  ■ . . ,  Wn/2  such  that 

\  '1/2 

/  w[t)  ■  (t)i{t)  dt  =  Y^  Wj  •  (t>^{tj),  (12) 

-1 

for  all  i  =  l,2,.'..,n,  with  (t)i,4>2,  ■  ■  ■  ,<i>n  the  first  n  eigenfunctions  of  the  operator  Q 
defined  in  (9). 

Corollary  2.3  The  above  theorem  provides  a  tool  for  the  efficient  approximate  evalua¬ 
tion  of  integrals  of  the  form  (12),  as  follows.  Given  a  positive  real  e,  we  construct  the 
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Singular  Value  Decomposition  of  the  operator  P  defined  in  (8).  Choosing  n  to  be  the 
smallest  even  integer  such  that 

OO 

E  (w) 

j=n+l 

we  construct  an  n/2-point  quadrature  that  integrates  n  first  right  singular  functions  ex¬ 
actly  ( effective  numerical  schemes  for  the  construction  of  such  quadratures  can  be  found 
in  [13,  28,  2]).  Now,  we  observe  that  due  to  the  triangle  inequality  combined  with  the 
positivity  of  the  obtained  weights  Wi,W2,  ■ . . ,  Wn/2, 
n/2  j 

\J2wj-  -  /  w{x)  ■  dt\  <  e  (14) 

j=i  -I 

for  any  X  G  [—1, 1]. 

Remark  2.2  The  principal  subject  of  this  note  is  the  fact  that  the  pattern  of  an  antenna 
array  is  formed  by  a  physical  process  amounting  to  a  hardware  implementation  of  a 
quadrature  formula  for  functions  of  the  form  (9).  Thus,  designing  a  configuration  of 
elements  for  such  an  antenna  is  equivalent  to  constructing  a  quadrature  formula  for 
functions  of  the  form(  9),  and  can  be  achieved  via  the  techniques  described  in  [13,  28,  2]). 

3  Elements  of  Antenna  Theory 

In  this  section,  we  summarize  certain  facts  about  the  theory  of  linear  antenna  arrays;  all 
of  these  facts  are  well-known,  and  can  be  found,  for  example,  in  [9]. 

3,1  Pattern  of  a  linear  array 

A  source  distribution  a  on  the  interval  [—1,1]  creates  the  far-field  pattern  /  :  [0,  tt]  — >•  (D 
given  by  the  formula 

1 

f{e)  =  J  a{u)  ■  du,  (15) 

-1 
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where  k  is  the  free-space  wavenumber,  u  is  the  point  on  the  interval  [—1,1],  and  6  is  the 
angle  between  the  point  on  the  horizon  where  the  far  field  is  being  evaluated  and  the 
rc-axis.  It  is  customary  to  introduce  the  notation 

X  =  cos{6),  (16) 

and  define  the  function  F  :  [— 1, 1]  — >  (D  by  the  formula 

F{x)  =  f{acos{x)).  (17) 

Now,  defining  the  operator  A  :  L‘^[—l,  1]  — >  1]  by  the  formula 

1 

A{a){x)  =  I  (j(u)  ■  du,  '  (18) 


we  observe  that 


1 

F  =  A{a)  =  I  a{u)  ■  du.  (19) 

-1 

The  function  F  is  usually  more  convenient  to  work  with  than  /,  and  the  following  obvious 
lemma  is  the  principal  reason  for  this  difference. 

Lemma  3.1  Suppose  that  a  G  L^[— 1,1],  the  function  F  G  L^]— 1,1]  is  defined  by  (19), 
a  is  a  real  number,  and  the  function  d  G  L^]— 1, 1]  is  defined  by  the  formula 

d{u)  =  •  a{u).  (20) 


Then 

A{a){x)  =  A{a){x  —  a)  (21) 

for  all  x  G  (— co,  oo).  In  other  words,  in  order  to  translate  the  antenna  pattern  F  (viewed 
as  a  function  of  x  =  cos{6)  )  by  a,  one  has  to  multiply  by  “  the  source  distribution  a 
generating  the  pattern  F. 


Observation  3.1  While  the  obvious  physical  considerations  lead  to  the  antenna  pattern 
F  defined  on  the  interval  [—1,1],  the  formulae  (15),  (17)  also  define  naturally  the  exten¬ 
sion  of  F  to  the  function  H  C;  in  a  mild  abuse  of  notation,  we  will  be  denoting  by  F 
both  the  original  mapping  [—1,1]  (D  and  its  extension  to  the  mapping  R,  — >  (D.  Simi¬ 
larly,  we  will  be  denoting  by  A  both  the  operator  L^[— 1,1]  — >•  L^[— 1,1]  defined  by  (18) 
and  its  natural  extension  mapping  L^[— 1, 1]  c°°{Il).  The  restriction  of  F  on  R\[— 1, 1] 

is  referred  to  as  the  invisible  spectrum  of  the  source  distribution  a  and  plays  an  important 
role  in  the  antenna  theory  (this  role  is  discussed  briefly  in  the  following  subsection).  By 
the  same  token,  the  restriction  of  F  on  the  interval  [—1,1]  is  referred  to  as  the  visible 
spectrum. 

When  an  antenna  array  is  implemented  in  hardware,  it  is  (usually)  constructed  of 
a  finite  collection  of  elements,  as  opposed  to  being  a  continuous  source  distribution. 
Mathematically,  it  is  equivalent  to  replacing  the  general  function  a  in  (15),  (19)  with  a 
defined  by  the  expression 

a{x)  =  (22) 

with  (fi,  (j)2, . . . ,  fin  the  source  distributions  generated  by  individual  elements,  and  the 
coefficients  fii,  P21  ■  ■  ■ ,  fln  the  intensities  of  the  elements.  As  a  rule,  the  elements  are 
localized  in  space  (i.e.  the  functions  fii,fi2,---,fin  are  supported  on  small  subintervals 
of  [—1, 1]),  and  very  often,  all  of  the  elements  are  identical  (i.e.  the  functions  fij  are 
translates  of  each  other),  so  that 

fij{u)  —  fi{u  —  Uj),  (23) 

with  fi  the  source  distribution  of  a  single  element  located  at  the  point  u  =  0,  and  Uj  the 
location  of  the  element  number  j.  Obviously,  the  far-field  pattern  of  fi  is  given  by  the 
formula 

1 

F^{x)  =  I  fi{u)  ■  du-  (24) 

-1 
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combining  (24)  with  (22)  and  (23),  we  obtain  the  identity 

)  n 

a{x)  =  /  (l>{u)  ■  du  (25) 

-1  i=i 

known  in  the  antenna  theory  as  the  principle  of  pattern  multiplication. 

Remark  3.2  The  standard  form  of  the  principle  of  multiplication  reads:  “The  field 
pattern  of  an  array  of  nonisotropic  but  similar  point  sources  is  the  product  of  the  pattern 
of  the  individual  source  and  the  the  pattern  of  an  array  of  isotropic  point  sources,  having 
the  same  locations,  relative  amplitudes  and  phases  as  the  nonisotropic  point  sources”  (see 
[9]).  Needless  to  say,  this  is  a  special  case  of  the  well-known  theorem  from  the  theory  of 
the  Fourier  Transform,  stating  that  the  Fourier  transform  of  the  product  of  two  functions 
is  the  convolution  of  the  Fourier  Transforms  of  multiplicants. 

4  Antenna  Patterns  and  Corresponding  Optimal  El¬ 
ement  Distributions 

4.1  Characteristics  of  an  antenna  pattern 

Depending  on  the  situation,  the  design  of  an  antenna  array  attempts  to  optimize  certain 
characteristics  of  the  resulting  far-held  pattern,  subject  to  certain  constraints  on  the 
number,  power,  etc.  of  the  elements.  Since  the  principal  purpose  of  this  note  is  to 
describe  a  technique  for  the  selection  of  the  locations  of  the  elements  that  approximate  a 
user-specihed  pattern,  we  could  use  any  reasonable  far-held  pattern  to  be  approximated. 
In  subsection  4.2,  4.3,  we  construct  optimal  element  distributions  for  the  so-called  sector 
patterns  and  cosecant  pattern,  respectively;  a  detailed  discussion  of  these  (and  several 
other)  pattern  cans  be  found,  for  example  in  [14]. 

We  will  say  that  the  antenna  pattern  has  the  e-bandwidth  b  if 

1 

J  \F{x)f  dx  =  ■  J  \F{x)\'^  dx  (26) 

6<||i||<l  —1 
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in  other  words,  the  proportion  of  the  energy  radiated  outside  the  e-beamwidth  from  the 
axis  of  the  beam  is  equal  to  e.  The  supergain  of  an  antenna  is  defined  (see,  for  example, 
[27]),  as  the  ratio 


4-00 

f  lF(a:)p  dx 

—  OO _ 

/  \F{xW  dx 
-1 


(27) 


The  supergain  (sometimes  referred  to  as  superdirectivity)  measures  the  ratio  of  the  en¬ 
ergy  associated  with  the  total  spectrum  of  the  antenna  to  the  energy  in  its  visible  spec¬ 
trum;  while  detailed  discussion  of  supergain  and  related  issues  is  outside  the  scope  of  this 
note,  we  will  observe  that  antenna  arrays  with  large  degrees  of  supergain  would  violate 
the  uncertainty  principle,  and  thus  are  physically  impossible.  Attempts  to  construct 
supergain  antennae  result  in  rapidly  (exponentially)  growing  Ohmic  losses,  prohibitive 
accuracy  requirements,  extremely  low  bandwidth,  etc.  Thus,  any  potentially  useful  pro¬ 
cedure  for  the  design  of  antenna  arrays  has  to  limit  the  supergain  of  the  resulting  patterns. 

4.2  Sector  patterns 

It  is  often  desirable  to  construct  antenna  patterns  that  are  as  constant  as  possible  within 
the  main  beam,  and  as  small  as  possible  outside  it;  in  other  words,  ideally,  the  pattern 
would  be  defined  by  the  formulae 


Fb{x)  =  1  for  |a;|  <  b, 
Fh{x)  =  0  for  |a;|  >  b, 


(28) 

(29) 


with  b  a  real  number  such  that  0  <  b  <  k.  Needless  to  say,  the  function  Fh  defined  by 
the  formulae  (28),  (29)  is  not  band-limited,  and  some  approximation  has  to  be  used.  A 
standard  procedure  is  to  truncate  the  Fourier  Transform  of  Fb,  approximating  it  by  the 
function  Fb  defined  by  the  formula 

sin{b  ■  t) 


>(^)  =  £ 


Akx-t 


(30) 
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(see,  for  example,  [26]).  An  important  special  case  occurs  when  b  =  k,  with  (30)  assuming 
the  form 


= £ 


sinjk  •  t)  _  i-k-x-t. 


(31) 


obviously,  the  latter  expression  is  a  band-limited  approximation  of  the  5-function.  An¬ 
other  frequently  encountered  situation  is  that  of  b  —  k/2,  so  that  (30)  assumes  the  form 


=  £ 


sm(|  •  t) 


Akxt 


(32) 


-1  t 

which  is  a  band-limited  approximation  to  the  beam  that  is  equal  to  1  for  —1/2  <  x  <  1/2 
and  to  zero  elsewhere. 

In  Section  4.4  below,  we  demonstrate  optimal  element  configurations  that  produce 
approximations  to  the  patterns  (31),  (32)  with  k  =  207r,  IOtt,  32.4676?:. 


Remark  4.1  While  (30)  is  by  no  means  the  only  possible  band-limited  approximations 
to  to  Fb^  it  is  quite  satisfactory  in  most  cases,  in  addition  to  being  simple.  Furthermore, 
the  principal  purpose  of  this  note  is  to  describe  a  technique  for  the  selection  of  locations 
of  the  nodes,  given  a  pattern  to  be  approximated.  Thus,  we  ignore  the  issue  of  the 
optimal  choice  of  Ft- 

4.3  Cosecant  patterns 

Another  standard  far-field  radiation  pattern  is  the  so-called  cosecant  pattern  (see,  for 
example,  [19]).  Given  two  real  numbers  0  <  a  <  6  <  1,  the  cosecant  pattern  Fa,b  is 
defined  by  the  formula 

Fa,bi^)  =  ^  (33) 

for  all  X  G  [a,  b],  and 

Fa,b{^)  =  0  (34) 
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for  all  X  e  ([—1, 1]  \  [fl,  6]).  Again,  the  function  Fa^b  defined  by  the  formulae  (33),  (34)  is 
not  band-limited,  and  can  not  be  represented  by  the  expression  of  the  form  (24).  Before 
the  scheme  of  this  note  can  be  applied  to  Fa^b,  the  latter  has  to  be  approximated  with  a 
band-limited  function;  as  discussed  in  Section  4.1  above,  if  such  an  approximation  is  to 
be  useful  as  an  antenna  pattern,  its  supergain  factor  has  to  be  controlled.  Fortunately, 
a  procedure  for  such  an  approximation  has  been  in  existence  for  more  than  35  years 
(see,  [18]);  the  algorithm  of  [18]  is  a  modification  of  the  least-squares  approach  permitting 
the  user  to  limit  the  supergain  factor  of  the  obtained  pattern  explicitly.  At  the  time,  the 
utility  of  the  scheme  of  [18]  was  limited  by  the  (perceived)  difficulty  in  the  numerical 
evaluation  of  Prolate  Spheroidal  Wave  functions;  given  the  present  state  of  numerical 
analysis,  this  difficulty  is  non-existent,  and  it  is  this  author’s  impression  that  the  insights 
of  [18],  [19]  deserve  more  attention  than  they  have  been  receiving. 

4.4  Optimal  distributions  of  elements 

In  this  subsection,  we  briefly  describe  an  algorithm  for  the  construction  of  optimal  (in 
the  sense  defined  below)  element  configurations  for  the  generation  of  antenna  patterns 
given  by  (15),  of  which  the  patterns  (29)-(31)  are  special  cases.  As  will  be  seen,  the 
procedure  is  in  fact  applicable  to  the  design  of  element  configurations  for  very  general 
far-field  patterns. 

We  start  with  observing  that  (15)  expresses  the  far-field  pattern  F  as  an  integral  over 
the  interval  [—1,1]  of  functions  of  the  form 

a{u)  ■  (35) 

with  X  =  cos{6)  determined  by  the  direction  9  in  which  the  far-field  is  being  evaluated.  In 
other  words,  the  problem  of  finding  efficient  antenna  element  distributions  is  equivalent 
to  that  of  constructing  quadrature  formulae  for  integrals  of  the  form  (8),  with 

w{t)  =  a{t).  (36) 
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In  the  cases  when  a  is  non-negative  everywhere  on  the  interval  [—1,1],  Theorem  2.2 
guarantees  the  existence  of  Generalized  Gaussian  Quadratures,  and  [13,  28])  provide  a 
satisfactory  numerical  apparatus  for  the  construction  of  such  quadratures.  Obviously,  the 
patterns  given  by  the  formula  (28)  are  not  generated  by  non-negative  source  distributions, 
except  when 

6  <  TT.  (37) 

Thus,  for  these  (and  many  other)  patterns,  the  conditions  of  Theorem  2.2  are  violated, 
and  the  existence  of  Generalized  Gaussian  Quadratures  is  not  guaranteed.  In  our  numer¬ 
ical  experiments,  the  techniques  of  [2])  (after  some  tuning)  have  always  been  successful 
in  finding  the  Gaussian  quadratures  for  integrals  of  the  form  (28);  some  of  our  results 
are  presented  in  Section  5  below. 

5  Numerical  Examples 

In  this  section,  we  present  examples  of  optimal  element  distributions  generating  the 
patterns  of  the  preceding  Section;  all  of  the  results  presented  here  have  been  obtained 
numericall3c  Antenna  patterns  we  present  are  compared  to  the  antenna  patterns  given 
by  uniform  source  distributions;  configurations  of  elements  approximating  these  antenna 
patterns  are  compared  to  equispaced  distributions  of  elements  generating  the  same  an¬ 
tenna  patterns. 

5.1  Optimal  distributions  of  elements 

In  this  section,  we  demonstrate  the  results  of  the  application  of  the  techniques  of  Sec¬ 
tion  4.4  of  this  note  to  the  types  of  antenna  patterns  described  in  the  Sections  4.2,  4.3. 

In  all  cases,  we  choose  the  size  of  an  antenna  array  and  a  pattern  to  be  reproduced,  and 
use  the  scheme  outlined  in  Section  4.4  to  design  a  distribution  of  antenna  elements  (both 
the  locations  and  the  intensities)  located  within  the  chosen  array  that  reproduces  the 
required  pattern.  For  comparison,  we  also  generate  optimal  (in  the  least  squares  sense) 


approximations  to  the  desired  pattern  generated  by  equispaced  elements  located  within 
the  same  array.  Since  the  number  of  equispaced  nodes  required  to  obtain  a  reasonable 
approximation  to  the  desired  pattern  is  (in  many  cases)  much  greater  than  the  number  of 
optimally  chosen  nodes,  for  each  example  we  demonstrate  patterns  generated  by  several 
such  configurations.  In  this  manner,  the  numbers  of  optimally  chosen  nodes  necessary 
to  obtain  reasonable  approximations  to  the  desired  patterns  can  be  compared  to  the 
numbers  of  equispaced  nodes  required  to  obtain  similar  results. 

5.1.1  Sector  patterns 

Example  5.1  The  first  example  we  consider  is  of  the  pattern  defined  by  the  formula  (32), 
with  k  =  62.8312,  so  that  the  size  of  the  array  is  20  wavelengths. 

In  Figure  5,  we  display  an  approximation  to  the  pattern  obtained  with  19  elements, 
overlayed  with  the  exact  pattern;  the  locations  of  the  elements  are  displayed  in  Figure  5a; 
the  relative  error  of  the  obtained  approximation  is  5.01%. 

Similarly,  in  Figure  5g,  we  display  the  approximation  to  the  pattern  obtained  with  21 
elements,  overlayed  with  the  exact  pattern;  the  relative  error  of  the  obtained  approxima¬ 
tion  is  0.443%;  in  Figure  5h,  we  display  the  the  approximation  obtained  with  17  elements. 
In  the  latter  case,  the  relative  error  of  the  obtained  approximation  is  6.43%;  Figure  5i 
depicts  the  17-node  distribution  producing  the  approximation  illustrated  in  Figure  5h. 
Finally,  Figure  5j  contains  a  graph  of  the  values  of  the  sources  located  at  the  17  nodes 
depicted  in  Figure  5i  and  generating  the  pattern  shown  in  Figure  5h. 

For  comparison,  the  optimal  approximation  obtained  with  19,  24,  29,  31,  and  34 
equispaced  elements  are  displayed  in  Figures  5b,  5c,  5d,  5e,  5f,  respectively;  these  are 
also  overlayed  with  the  exact  pattern. 

Example  5.2  Our  second  example  is  identical  to  the  first  one,  with  the  exception  that 
k  =  31.416,  so  that  the  size  of  the  array  is  10  wavelengths. 

In  Figure  6,  we  display  an  approximation  to  the  pattern  obtained  with  9  elements, 
overlayed  with  the  exact  pattern;  the  locations  of  the  elements  are  displayed  in  Figure  6a; 
the  relative  error  of  the  obtained  approximation  is  11.2%. 


Similarly,  in  Figure  6f,  we  display  the  approximation  to  the  pattern  obtained  with  11 
elements,  overlayed  with  the  exact  pattern;  the  relative  error  of  the  obtained  approxima¬ 
tion  is  0.600%. 

For  comparison,  the  optimal  approximation  obtained  with  9,  14,  16,  and  18  equispaced 
elements  are  displayed  in  Figures  6b,  6c,  6d,  5e,  respectively;  these  are  also  overlayed 
with  the  exact  pattern. 

Example  5.3  Our  third  example  is  identical  to  the  preceding  two,  with  the  exception 
that  k  =  102,  so  that  the  size  of  the  array  is  about  32.45  wavelengths. 

In  Figure  7a,  we  display  an  approximation  to  the  pattern  obtained  with  23  optimally 
distributed  elements,  overlayed  with  the  exact  pattern  and  with  the  pattern  obtained  with 
23  equispaced  elements. 

The  relative  error  of  the  obtained  approximation  is  5.4%;  needless  to  say,  the  error  of 
the  approximation  obtained  with  the  equispaced  nodes  is  more  than  70%.  As  can  be  seen 
from  Figure  7c,  the  actual  size  of  the  obtained  2Z-element  array  is  about  21  wavelengths; 
in  other  words,  in  order  to  obtain  this  precision,  the  array  needs  to  be  about  2/5  of  the 
nominal  (maximum  permitted)  length. 

In  Figure  7b,  we  display  the  approximation  to  the  pattern  obtained  with  42  and  48 
elements,  overlayed  with  the  exact  pattern. 

It  is  worth  noting  that  with  33  optimally  distributed  elements,  the  pattern  is  approxi¬ 
mated  to  the  precision  0.12%;  we  do  not  display  the  obtained  pattern  since  it  is  visually 
indistinguishable  from  the  pattern  being  approximated. 

Example  5.4  Our  final  example  is  somewhat  different  from  the  preceding  ones,  in  that 
instead  of  approximating  a  sector  pattern,  we  approximate  a  cosecant  pattern  (see  (33),  (3/) 
in  Subsection  f.3  above). 

In  this  example,  we  set 

a  =  sm(15°),  (38) 
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b  =  sin{75°), 


(39) 


and  use  the  procedure  of  [18]  to  approximate  Fa,b  with  a  band-limited  function.  The  band- 
limit  has  been  more  or  less  arbitrarily  set  to  110,  resulting  in  an  antenna  array  about  35 
wavelengths  in  size,  and  the  supergain  factor  of  the  approximation  was  set  to  1.1. 

In  Figure  8a,  we  display  an  approximation  to  the  pattern  obtained  with  53  optimally 
distributed  elements,  overlayed  with  the  exact  bandlimited  pattern  and  with  the  pattern 
obtained  with  53  equispaced  elements. 

The  relative  error  of  the  obtained  approximation  is  1.79%;  the  error  of  the  approxi¬ 
mation  obtained  with  the  equispaced  nodes  is  about  42%. 

In  Figure  8b,  we  display  the  approximation  to  the  pattern  obtained  with  47  optimally 
distributed  elements,  overlayed  with  the  exact  pattern;  the  purpose  of  this  final  figure  is 
to  demonstrate  the  behavior  of  the  scheme  when  the  number  of  elements  is  insufficient 
(i.e.  when  the  array  is  underresolved). 

It  is  worth  noting  that  it  takes  about  70  equispaced  nodes  to  obtain  the  resolution 
obtained  with  47  optimally  chosen  ones. 

The  following  observations  can  be  made  from  Figures  5  -  8b,  and  from  the  more 
detailed  numerical  experiments  performed  by  the  author. 

1.  In  order  to  obtain  reasonable  precision,  the  scheme  requires  about  1  point  per  wave¬ 
length  in  the  antenna  array;  this  is  more  or  less  independent  from  the  structure  of  the 
beam  as  long  as  the  pattern  is  symmetric  about  the  point  a:  =  0.  This  fact  is  observed 
numerically,  even  for  modest  numbers  of  nodes;  for  large-scale  arrays,  this  statement 
(interpreted  asymptotically)  can  be  proved  rigorously.  For  certain  beam  structures,  the 
required  number  of  nodes  is  even  less  (see  Example  5.3).  The  reasons  for  these  additional 
savings  are  subtle,  and  have  to  do  with  the  fact  that  the  continuous  source  distribution 
generating  the  pattern  is  relatively  small  on  a  large  part  of  the  antenna  array;  the  al¬ 
gorithm  of  [2]  takes  advantage  of  this  fact  to  reduce  the  number  of  nodes.  When  the 
beam  is  not  symmetric  about  a:  =  0,  the  number  of  elements  required  does  depend  on 
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the  structure  of  the  pattern,  and  the  dependence  is  fairly  complicated.  Generally,  the 
improvement  for  non-symmetric  beams  is  less  than  that  for  the  symmetric  ones. 

2.  The  qualiative  behavior  of  the  scheme  is  similar  to  that  of  the  Gaussian  quadratures 
in  that  it  displays  no  convergence  at  all  until  a  certain  minimum  number  of  nodes  is 
achieved;  after  that,  the  convergence  is  very  fast.  This  behavior  is  not  surprising,  since 
the  scheme  is  based  on  a  Generalized  Gaussian  quadrature. 

3.  For  the  sector  pattern  with  the  sector  [—1/2, 1/2],  the  scheme  reduces  the  required 
number  of  nodes  by’ a  factor  of  about  1.5  for  small-scale  problems,  and  roughly  by  a 
factor  of  2  for  large-scale  ones;  again,  for  large-scale  problems,  an  asymptotic  version  of 
this  statement  can  be  proven  rigorously. 

4.  For  the  cosecant  pattern  with  the  parameters  specified  by  (38),  (39),  the  number 
of  nodes  required  is  reduced  by  approximately  a  factor  of  1.4.  As  the  sidelobe  level  is 
reduced,  the  improvement  obtained  by  going  from  the  equispaced  discretization  to  the 
optimal  one  increases  rapidly. 

5.  An  examination  of  Figures  5a,  6a  shows  that  w'hile  the  optimal  nodes  are  by  no  means 
uniform,  they  display  no  clustering  behavior. 

6.  An  examination  of  Figure  5j  shows  that  the  intensities  of  individual  elements  do  not 
become  large;  this  is  confirmed  by  the  more  extensive  numerical  experiments  performed 
by  the  author. 

7.  The  combination  of  the  preceding  two  paragraphs  (combined  with  additional  numer¬ 
ical  experiments  and  analysis)  provide  evidence  that  configurations  of  this  type  should 
pose  no  supergain  problems. 

6  Generalizations 

The  results  described  above  admit  radical  generalizations  in  several  directions;  several 
such  directions  are  discussed  below, 


1.  Conformal  one-dimensional  arrays.  The  extension  of  the  techniques  of  this  note 
to  one-dimensional  arrays  located  on  curves  in  is  completely  straightforward,  involving 
only  a  modest  increase  of  the  CPU  time  requirements  of  the  procedure.  Improvement  in 
the  number  of  nodes  required  to  produce  a  prescribed  pattern  is  similar  to  that  in  the 
case  of  a  linear  array. 

2.  Planar  two-dimensional  arrays.  A  straightforward  generalization  of  the  results  of 
Sections  4,  5,  is  to  rectangular  planar  arrays.  Here,  a  tensor  product  quadrature  can  be 
constructed  from  the  quadratures  of  Sections  4,  5,  possessing  all  of  the  desirable  prop¬ 
erties  of  the  latter.  Obviously,  the  advantage  in  the  number  of  transducers  is  squared, 
so  that  (for  example)  replacing  50  nodes  in  each  of  the  two  directions  by  23  nodes  (see 
Example  5.3  above)  will  lead  to  a  factor  of  (50/23)^  ~  4.7  savings  in  the  number  of 
elements. 

The  theory  of  Section  4  has  been  extended  for  disk-shaped  arrays,  via  {inter  alia)  the 
techniques  developed  in  [23].  The  improvement  in  the  number  of  nodes  is  comparable  to 
that  obtained  in  the  rectangular  geometry,  and  the  CPU  time  requirements  do  not  differ 
appreciably  from  those  in  the  case  of  linear  one-dimensional  arrays. 

The  extension  of  the  theory  to  more  general  geometries  in  the  plane  is  in  progress.  At 
the  present  time,  our  only  numerical  experiments  have  been  with  arrays  on  triangles;  the 
results  are  encouraging,  but  the  CPU  time  requirements  of  the  algorithms  are  excessive 
(we  have  only  been  able  to  design  triangular  arrays  about  6  wavelengths  in  size).  We 
are  now  in  the  process  of  constructing  a  more  efficient  numerical  procedure  for  such 
computations. 

3.  Conformal  two-dimensional  arrays.  The  only  environment  in  which  we  have 
a  satisfactory  theory  is  when  the  array  is  located  on  a  surface  of  revolution;  even  in 
this  environment,  no  experiments  have  been  performed.  We  have  not  investigated  more 
general  conformal  two-dimensional  arrays  in  sufficient  detail. 
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Figure  5:  The  pattern  created  by  the  19  optimal  elements,  depicted  in  Figure 
5a  as  described  in  Example  5.1 
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Figure  5a:  The  distribution  of  elements  creating  the  pattern  depicted  in 
Figure  5,  as  described  in  Example  5.1 
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Figure  5c:  The  optimal  approximation  to  the  sector  pattern  generated  by  24 
equispaced  nodes,  as  described  in  Example  5.1 
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Figure  5d:  The  optimal  approximation  to  the  sector  pattern  generated  by  29 
equispaced  nodes,  as  described  in  Example  5.1 


Figure  5e:  The  optimal  approximation  to  the  sector  pattern  generated  by  31 
equispaced  nodes,  as  described  in  Example  5.1 
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Figure  5f:  The  optimal  approximation  to  the  sector  pattern  generated  by  34 
equispaced  nodes,  as  described  in  Example  5.1 


Figure  5g:  The  optimal  approximation  to  the  sector  pattern  generated  by  21 
optimal  nodes,  as  described  in  Example  5.1 
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Figure  6b:  The  optimal  approximation  to  the  sector  pattern  generated  by  9 
equispaced  nodes,  as  described  in  Example  5.2 


Figure  6c:  The  optimal  approximation  to  the  sector  pattern  generated  by  14 
equispaced  nodes,  as  described  in  Example  5.2 
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Figure  6d:  The  optimal  approximation  to  the  sector  pattern  generated  by  16 
equispaced  nodes,  as  described  in  Example  5.2 


Figure  6e:  The  optimal  approximation  to  the  sector  pattern  generated  by  18 
equispaced  nodes,  cis  described  in  Example  5.2 


Figure  6f:  The  pattern  created  by  the  11  optimal  elements,  in  Example  5.2 
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Figure  7a:  The  approximation  to  the  sector  pattern  generated  by  23  optimal 
elements,  vs.  optimal  approximation  by  23  equispaced  nodes,  as  described  in 

Example  5.3 
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Figure  8a:  The  approximation  to  the  cosecant  pattern  generated  by  53 
optimal  elements,  vs.  optimal  approximation  by  53  equispaced  nodes,  as 
described  in  Example  5.4 
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Figure  8a:  The  approximation  to  the  cosecant  pattern  generated  by  47 
optimal  elements,  as  described  in  Example  5.4 
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Polynomials  are  one  of  principal  tools  of  classical  numerical  analysis.  When  a  function 
needs  to  be  interpolated,  integrated,  differentiated,  etc.,  it  is  assumed  to  be  approximated 
by  a  polynomial  of  a  certain  fixed  order  (though  the  polynomial  is  almost  never  constructed 
explicitly),  and  a  treatment  appropriate  to  such  a  polynomial  is  applied.  We  introduce  anal¬ 
ogous  techniques  based  on  the  assumption  that  the  function  to  be  dealt  with  is  band-limited, 
and  use  the  well-developed  apparatus  of  Prolate  Spheroidal  Wave  Functions  to  construct 
quadratures,  interpolation  and  differentiation  formulae,  etc.  for  band-limited  functions. 
Since  band-limited  functions  are  often  encountered  in  physics,  engineering,  statistics,  etc. 
the  apparatus  we  introduce  appears  to  be  natmal  in  many  environments.  Our  results  are 
illustrated  with  several  numerical  examples. 
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1  Introduction 

Numerical  quadrature  and  interpolation  are  a  well-developed  part  of  numerical  analysis; 
polynomials  are  the  classical  tool  for  the  design  of  such  schemes.  Conceptually  speaking, 
one  assumes  that  the  function  is  well-approximated  by  expressions  of  the  form 

(1) 

j=0 

with  reasonably  small  n,  and  designs  algorithms  that  are  effective  for  functions  of  the 
form  (1)  (needless  to  say,  one  almost  never  actually  computes  the  coefficients  {oj};  one 
only  uses  the  fact  of  their  existence).  Obviously,  the  polynomial  approach  is  only  effective 
for  functions  that  are  well- approximated  by  polynomials. 

When  one  has  to  handle  functions  that  are  well-behaved  on  the  whole  line  (for  ex¬ 
ample,  in  signal  processing),  polynomials  are  not  an  appropriate  tool.  In  such  cases, 
trigonometric  polynomials  are  used;  existing  tools  are  very  satisfactory  for  dealing  with 
functions  defined  and  well-behaved  on  the  whole  of  Such  tools,  in  effect,  make  the 
assumption  that  the  functions  are  band-limited  or  nearly  so;  a  function  /  :  R  — >■  R  is 
said  to  be  band-limited  if  there  exist  a  positive  real  c  and  a  function  a  6  L^[-l,  1]  such 
that 

fix)  =  |_%'“V(t)dt.  (2) 

However,  in  many  cases,  we  are  confronted  with  band-limited  functions  defined  on  inter¬ 
vals  (or,  more  generally,  on  compact  regions  in  R”).  Wave  phenomena  are  a  rich  source 
of  such  functions,  both  in  the  engineering  and  computational  contexts;  they  are  also 
encountered  in  fluid  dynamics,  signal  processing,  and  many  other  areas.  Often,  such 
functions  can  be  effectively  approximated  by  polynomials  via  standard  tools  of  classical 
analysis.  However,  even  when  such  approximations  are  feasible,  they  are  usually  not 
optimal.  Smooth  periodic  functions  are  a  good  illustration  of  this  observation:  while 
they  can  be  approximated  by  polynomials  (for  example,  via  Chebyshev  or  Legendre 
expansions),  they  are  more  efficiently  approximated  by  Fourier  expansions,  both  for  an¬ 
alytical  and  numerical  purposes.  It  would  appear  that  an  approach  explicitly  based  on 
trigonometric  polynomials  could  be  more  efficient  in  dealing  with  band-limited  functions. 

In  the  engineering  context,  such  an  apparatus  was  constructed  more  than  30  years 
ago  (see  [20]-[21],  [7]-[9]).  The  natural  tool  for  analyzing  band-limited  functions  on  R^  is 
the  Fourier  Transform,  unless  the  functions  are  periodic,  in  which  case  the  natural  tool  is 


the  Fourier  Series.  The  authors  of  [20]-[21]  observe  that  for  the  analysis  of  band-limited 
functions  on  the  interval,  Prolate  Spheroidal  Wave  Functions  are  likewise  a  natural  ap¬ 
proach.  The  authors  also  construct  a  multidimensional  version  of  the  theory,  though 
their  apparatus  is  only  complete  for  the  case  of  spherical  regions. 

The  present  paper  constructs  tools  for  the  use  of  the  approach  of  [20]- [21]  in  the 
modern  computational  environment.  We  construct  a  class  of  quadratures  for  band- 
limited  functions  that  closely  parallel  the  Gaussian  quadratures  for  polynomials.  The 
nodes  are  very  close  to  being  roots  of  appropriately  chosen  Prolate  Spheroidal  Wave 
Functions,  the  resulting  quadratures  are  stable,  and  all  weights  are  positive.  As  in  the 
case  of  polynomials,  there  are  interpolation,  differentiation  and  indefinite  integration 
schemes  associated  with  the  obtained  quadratures,  exact  on  certain  classes  of  band- 
limited  functions.  These  procedures  are  the  main  tools  necessary  for  the  numerical  use 
of  spectral  discretizations  based  on  Prolate  Spheroidal  Wave  Functions,  instead  of  on 
the  usual  polynomial  bases.  When  dealing  with  band-limited  functions,  the  number  of 
nodes  required  by  these  procedures  to  obtain  a  prescribed  accuracy  is  much  less  than 
that  required  by  their  polynomial-based  counterparts.  An  additional  bonus  is  the  fact 
that  the  condition  number  of  differentiation  of  prolate  spheroidal  wave  functions  is  less 
than  that  of  differentiation  of  the  usual  polynomial  basis  functions  (see  Section  8  below). 

This  paper  is  organized  as  follows.  Section  2  summarizes  various  standard  mathemat¬ 
ical  facts  used  in  the  remainder  of  the  paper.  Section  3  contains  derivations  of  various 
results  used  in  the  algorithms  described  in  later  sections.  Section  4  describes  algorithms 
for  evaluation  of  prolate  spheroidal  wave  functions  and  associated  eigenvalues.  Section  5 
describes  a  construction  of  quadratures  for  band-limited  functions.  Section  6  describes 
an  alternative  approach  to  arriving  at  such  quadratures;  it  shows  that  roots  of  appropri¬ 
ately  chosen  prolate  spheroidal  wave  functions  can  serve  as  quadrature  nodes.  Section  7 
analyzes  the  use  of  prolate  spheroidal  wave  functions  for  interpolation.  Section  8  con¬ 
tains  results  of  our  numerical  experiments  with  quadratures  and  interpolation.  Section  9 
contains  a  number  of  miscellaneous  properties  of  prolate  spheroidal  wave  functions,  and 
Section  10  contains  generalizations  and  conclusions. 

2  Mathematical  Preliminaries 

As  a  matter  of  convention,  in  this  paper  the  norm  of  a  function  is,  unless  stated  otherwise, 
its  norm: 


ll/ll  =  \  dx. 
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2.1  Chebyshev  systems 

Definition  2.1  4  sequence  of  functions  4>\,---i4>n  will  be  referred  to  as  a  Chebyshev 
system  on  the  interval  [a,  6]  if  each  of  them  is  continuous  and  the  determinant 


M^i)  •••  (f>iM 

is  nonzero  for  any  sequence  of  points  Xi,...,Xn  such  that  a  <  Xi  <  X2 . . .  <  Xn  <  b. 


(4) 


An  alternate  definition  of  a  Chebyshev  system  is  that  any  linear  combination  of  the 
functions  with  nonzero  coefficients  must  have  fewer  than  n  zeros. 

Examples  of  Chebyshev  and  extended  Chebyshev  systems  include  the  following  (ad¬ 
ditional  examples  can  be  found  in  [11]). 


Example  2.1  The  powers  l,a:,a:^, . . . form  an  extended  Chebyshev  system  on  the 
interval  {—oo,oo). 

Example  2.2  The  exponentials  . . . ,  form  an  extended  Chebyshev  sys¬ 
tem  for  any  Ai, . . . ,  >0  on  the  interval  [0,  oo). 


Example  2.3  The  functions  1,  cos  x,  sin  x,  cos  2x,  sin  2a:, ... ,  cos  nx,  sin  nx  form  a  Cheby¬ 
shev  system  on  the  interval  [0,27r]. 


2.2  Generalized  Gaussian  quadratures 

A  quadrature  rule  is  an  expression  of  the  form 

(5) 

where  the  points  Xj  €  R,  and  coefficients  Wj  €  R,  are  referred  to  as  the  nodes  and  weights 
of  the  quadrature,  respectively.  They  serve  as  approximations  to  integrals  of  the  form 

rb 

/  ^{x)u{x)  dx,  (6) 

Ja 

with  oj  being  an  integrable  non-negative  function. 

Quadratures  are  typically  chosen  so  that  the  quadrature  (5)  is  equal  to  the  desired 
integral  (6)  for  some  set  of  functions,  commonly  polynomials  of  some  fixed  order.  Of 
these,  the  classical  Gaussian  quadrature  rules  consist  of  n  nodes  and  integrate  polynomi¬ 
als  of  order  2n  —  1  exactly.  In  [13],  the  notion  of  a  Gaussian  quadrature  was  generalized 
as  follows: 
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Definition  2.2  A  quadrature  formula  will  be  referred  to  as  Gaussi&n  with  respect  to  a 
set  of  2n  functions  ,  02n  :  [o-,  6]  — >  R.  and  a  weight  function  lo  :  [a,  b]  — )•  if  it 

consists  of  n  weights  and' nodes,  and  integrates  the  functions  4>i  exactly  with  the  weight 
function  u  for  all  i  =  1, . . . ,  2n.  The  weights  and  nodes  of  a  Gaussian  quadrature  will  be 
referred  to  as  Gaussian  weights  and  nodes  respectively. 

The  following  theorem  appears  to  be  due  to  Markov  [14,  15];  proofs  of  it  can  also  be 
found  in  [12]  and  [11]  (in  a  somewhat  different  form). 

Theorem  2.1  Suppose  that  the  functions  <f>i, . . . ,  (j>2n  •  [oj^]  R  /orm  a  Chebyshev 
system  on  [a,  6] .  Suppose  in  addition  that  a;  :  [a,  6]  — ^  R  is  a  non-negative  integrable 
function  [a,  6]  R.  Then  there  exists  a  unique  Gaussian  quadrature  for  the  functions 
^1)  •  •  • )  <f>2n  on  [a,  6]  with  respect  to  the  weight  function  u.  The  weights  of  this  quadrature 
are  positive. 

While  the  existence  of  Generalized  Gaussian  Quadratures  was  observed  more  than 
100  years  ago,  the  constructions  found  in  [14,  15],  [6,  12],  [10,  11]  do  not  easily  yield 
numerical  algorithms  for  the  design  of  such  quadrature  formulae;  such  algorithms  have 
been  constructed  recently  (see  [13,  25,  2]). 

Remark  2.1  It  might  be  worthwhile  to  observe  here  that  when  a  Generalized  Gaussian 
quadrature  is  to  be  constructed,  the  determination  of  its  nodes  tends  to  be  the  critical 
step  (though  the  procedure  of  [13,  25,  2]  determines  the  nodes  and  weights  simultane¬ 
ously).  Indeed,  once  the  nodes  xi,X2,...,Xn  have  been  found,  the  weights  Wi,W2,...,Wn 
can  be  determined  easily  as  the  solution  of  the  n  x  n  system  of  linear  equations 

•  0i(a;j)  =  /  (l)i{x)  dx,  (7) 

j=i 

with  i  =  1,  2, . . .  ,  n. 

2.3  Legendre  Polynomials 

In  agreement  with  standard  practice,  we  will  be  denoting  by  the  classical  Legendre 
polynomials,  defined  by  the  three-term  recursion 

Pn+l{x)  =  ^^TPn(x)--^Pr.-,{x),  (8) 

n  -h  1  n  +  1 

with  the  initial  conditions 

Po{x)  =  1,  (9) 

Pi(a:)  =  x] 
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as  is  well-known 


PkU)  =  1  (10) 

for  all  /:  =  0, 1, 2, . . and  each  of  the  polynomials  Pk  satisfies  the  differential  equation 

+  (11) 

The  polynomials  defined  by  the  formulae  (8), (9)  are  orthogonal  on  the  interval  [—1, 1]; 
however,  they  are  not  orthonormal,  since  for  each  n  >  0, 

=  ^2- 

the  normalized  version  of  the  Legendre  polynomials  will  be  denoted  by  P„,  so  that 

^(z)=Pn(rr)-yn  +  l/2.  (13) 

The  following  lemma  follows  immediately  from  the  Cauchy-Schwartz  inequality  and  from 
the  orthogonality  of  the  Legendre  polynomials  on  the  interval  [—1,1]: 

Lemma  2.2  For  all  integer  k>n, 


Pn{x)  dx 


< 


For  all  integer  0  <  k  <  n, 


(14) 


2.4  Convolutional  Volterra  Equations 

A  convolutional  Volterra  equation  of  the  second  kind  is  an  expression  of  the  form 

(p{x)  =  f  K (x  —  t)  (p{t)  dt  -i-  a{x)  (16) 

Ja 

where  a,  b  are  a  pair  of  numbers  such  that  a  <  b,  the  functions  a,  K  :  [a,  b]  are 
square-integrable,  and  (/?  ;  [a,  6]  — >•  C  is  the  function  to  be  determined.  Proofs  of  the 
following  theorem  can  be  found  in  [4],  as  well  as  in  many  other  sources. 

Theorem  2.3  The  equation  (16)  always  has  a  unique  solution  on  the  interval  [a, 6].  If 
both  functions  K,  o  are  k  times  continuously  differentiable,  the  solution  (p  is  also  k  times 
continuously  differentiable. 


2.5  Prolate  Spheroidal  Wave  Functions 

In  this  subsection,  we  summarize  certain  facts  about  the  Prolate  Spheroidal  Wave  Func¬ 
tions.  Unless  stated  otherwise,  all  these  facts  can  be  found  in  [20,  17]. 

Given  a  real  c  >  0,  we  will  denote  by  Fc  the  operator  X^[-l,  1]  ->■  X^[-l,  1]  defined 
by  the  formula 

F,(<p){x)  ^  ^(t)  dt.  (17) 

Obviously,  Fc  is  compact;  we  will  denote  by  Aq,  Aj, . . . ,  A„, . . .  the  eigenvalues  of  Fc 
ordered  so  that  |Aj_i|  >  [Ajl  for  all  natural  j.  For  each  non-negative  integer  j,  we  will 
denote  by  ipj  the  eigenfunctions  corresponding  to  Xj,  so  that 

Xjtjjj{x)  =  J  dt,  (18) 

for  all  rr  G  [—1,1];  we  adopt  the  convention  that  the  functions  are  normalized  such  that 
=  1,  for  all  j.^  The  following  theorem  is  a  combination  of  several  lemmas 
from  [20], [6], [11]. 

Theorem  2.4  For  any  positive  real  c,  the  eigenfunctions  ‘ipo,'ipi, ... ,  of  the  operator  Fc 
are  purely  real,  are  orthonormal,  and  are  complete  in  X^[— 1,1].  The  even-numbered 
eigenfunctions  are  even,  and  the  odd-numbered  ones  are  odd.  All  eigenvalues  of  Fc 
are  non-zero  and  simple;  the  even-numbered  eigenvalues  are  purely  real,  and  the  odd- 
numbered  ones  are  purely  imaginary;  in  particular,  Xj  =  i^Xj\.  The  functions  ipi  consti¬ 
tute  a  Chebychev  system  on  the  interval  [—1, 1];  in  particular,  the  function  ipi  has  exactly 
i  zeroes  on  that  interval,  for  any  i  =  0, 1, . . . ,. 

We  will  define  the  self-adjoint  operator  Qc :  X^[— 1, 1]  — )■  X^[— 1, 1]  by  the  formula 

^  /  N  1  7^  sinic  ■ix  —  t)),.,  , 

<3c(^)  =  -  /  - - — (19) 

TT  J-l  X  —  t 
a  simple  calculation  shows  that 

<3.  =  ■  i?  •  fc,  (20) 

that  Qc  has  the  same  eigenfunctions  as  Fc,  and  that  the  j-th  (in  descending  order) 
eigenvalue  p,j  of  Qc  is  connected  with  Xj  by  the  formula 

H,  =  ^-\Xj\\  (21) 

^This  convention  differs  from  that  used  in  [20];  however,  the  present  paper  is  concerned  almost 
exclusively  with  approximation  of  functions  on  [—1,1],  and  in  that  context,  the  convention  that  the 
functions  {tpj}  have  unit  norm  on  that  interval  is  by  far  the  most  convenient. 
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The  operator  Qc  is  obviously  closely  related  to  the  operator  Pc :  L^[— oo,  oo]  — >  [—00, 00] 
defined  by  the  formula 

„  ,  .  1  r°°  sin(c  ■  (x  —  t))  ,  .  , 

=  - T - -■(p{t)dt,  (22) 

TT  •/— 00  X  t 

which,  as  is  well  known,  is  the  orthogonal  projection  operator  onto  the  space  of  functions 
of  band  limit  c  on  (—00, 00). 

For  large  c,  the  spectrum  of  Qc  consists  of  three  parts:  about  2c/7c  eigenvalues  that 
are  very  close  to  1,  followed  by  order  log(c)  eigenvalues  which  decay  exponentially  from  1 
to  nearly  0;  the  remaining  eigenvalues  are  all  very  close  to  zero.  The  following  theorem, 
proven  (in  a  slightly  different  form)  in  [19],  describes  the  spectrum  of  Qc  more  precisely. 

Theorem  2.5  For  any  positive  real  c  and  0  <  a  <  1  the  number  N  of  eigenvalues  of 
the  operator  Qc  that  are  greater  than  a  satisfies  the  inequality 

—  +  log  - — log(c)  -  10  •  log(c)  <N  <  (23) 

TT  VTT'^  a  J 

—  +  log  ^ log(c)  +  10  •  log(c) . 

TT  VTT-^  a  J 

By  a  remarkable  coincidence,  the  eigenfunctions  ,ipn  of  the  operator  Qc  turn 

out  to  be  the  Prolate  Spheroidal  Wave  functions,  well-known  from  classical  Mathematical 
Physics  (see,  for  example,  [16]).  The  following  theorem  formalizes  this  statement;  it  is 
proven  in  a  considerably  more  general  form  in  [21]. 

Theorem  2.6  For  any  c>  0,  there  exists  a  strictly  increasing  sequence  of  positive  real 
numbers  X0)Xi)  •  •  •  such  that  for  each  j  >  0,  the  differential  equation 

(1  -  x'^)  ip"{x)  —  2x'tf'{x)  +  {xj  -  ^^)  =  0  (24) 

has  a  solution  that  is  continuous  on  the  interval  [—1, 1].  For  each  j  >  0,  the  function  'ipj 
(defined  in  Theorem  2.4)  is  the  solution  of  (24)- 

3  Analytical  Apparatus 

3.1  Prolate  Series 

Since  the  functions  '00;  '0i,  ■  •  - ,  -  are  a  complete  orthonormal  basis  in  L^[— 1, 1],  any 

formula  for  the  inner  product  of  prolate  spheroidal  wave  functions  with  another  function 
/  is  also  a  formula  for  the  coefficients  of  an  expansion  of  /  into  prolate  spheroidal  func¬ 
tions  (which  we  will  refer  to  as  the  prolate  expansion  of  /).  Thus  the  following  theorem 
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provides  the  coefficients  of  the  prolate  expansion  of  the  derivative  of  a  prolate  spheroidal 
function,  and  also  the  coefficients  of  the  prolate  expansion  of  a  prolate  spheroidal  wave 
function  multiplied  by  x.  Those  coefficients  are  also  the  entries  of  the  matrix  for  differen¬ 
tiation  of  a  prolate  expansion  (producing  another  prolate  expansion),  and  the  entries  of 
the  matrix  for  multiplication  of  a  prolate  expansion  by  2:,  respectively.  (These  formulae 
are  not,  however,  suitable  for  producing  such  matrices  numerically,  since  in  many  cases 
they  exhibit  catastrophic  cancellation.) 

Theorem  3.1  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  If  m  =  n  {mod  2) ,  then 

^Prnix)  dx  =  j  ^x  ^„(rr)  i)m{^)  dx  =  0.  (25) 

Ifm^n  (mod  2),  then 

j^^^'n{x)'(frnix)  dx  = 

J  ^X1pn{x)lpmix)  dx  = 

Proof.  Since  the  functions  'ifj  are  alternately  even  and  odd,  (25)  is  obvious.  In  order  to 
prove  (26),  we  start  with  the  identity 

An'^^n  =  i^nit)  dt  (28) 

(see  (18)  in  Subsection  2.5).  Differentiating  (28)  with  respect  to  x,  we  obtain 


2A2 


>■1  +  AS 

2  AtoA„ 


^m(l)l/’n(l), 


ic  A^  -h  A2 


^m(l)^n(l)- 


(26) 

(27) 


Xn1pn{x)  J  te*“‘V'n(i)  dt. 


(29) 


Projecting  both  sides  of  (29)  on  iprn  and  using  the  identity  (28)  (with  n  replaced  with 
m)  again,  we  have 


Xn  ll^mix)  dx 

=  ic  J  'ipm{x)  J  te*“*^„(t)  dtdx 

=  icj  tlpn{t)  J  ll’m{x)  dxdt 

—  i  C  Xjjri  J"  ^tlpji{t^'4^m{t)  dt. 


(30) 
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Obviously,  the  above  calculation  can  be  repeated  with  m  and  n  exchanged,  yielding  the 
identity 

combining  (30)  with  (31),  we  have 

dx  =  ^  j^^'tprn{x)'lpn{x)  dx.  (32) 

On  the  other  hand,  integrating  the  left  side  of  (32)  by  parts,  we  have 

^n{x)  dx 

=  -0771(1)  V'n(l)  -  '0^(-l)t/)„(-l)  -  j^^'i})'^{x)il)^{x)  dx.  (33) 

Since  m^n  (mod  2),  we  rewrite  (33)  as 
j  ^'ip'^ix)  iJn{x)  dx 

=  2xprn{l)i’n{^)  -  j_^'>P'n{x)'ipm{x)  dx.  (34) 

Now,  combining  (32)  and  (34)  and  rearranging  terms,  we  get 

fi  2 

/  i^m{x)  dx  =  3  ^m(l)  07i(l)  •  (35) 

•'-1 

Substituting  (30)  into  (35),  we  get 


j^^X'll;n{x)^m{x)  dx 
^  ^  /  xp'„{x) 'tpmix)  dx 

Atti  J 


1_K 

ic  Xm  +  A2 

2  XffiXji 
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Ml)  Ml) 


ic  A^  +  A2 


0’7n(l)07i(l)- 


(36) 

□ 


The  following  corollary,  which  is  an  immediate  consequence  of  (32),  finds  use  in  the 
numerical  evaluation  of  the  eigenvalues  {Aj}: 


Corollary  3.2  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  If m^n  (mod 2),  then 


A2  /  dx 

'^m  J  —  1 _ 


(37) 


3.2  Decay  of  Legendre  Coefficients  of  Prolate  Spheroidal  Wave- 
functions 

Since  each  of  the  functions  xpj  is  analytic  on  C,  on  the  interval  [—1,1]  it  can  be  expanded 
in  a  Legendre  series  of  the  form 

OO 

=  S  ^kPk{x),  (38) 

k=Q 

with  the  coefficients  /3fc  decaying  superalgebraically;  the  following  two  theorems  establish 
bounds  for  the  decay  rate. 


Lemma  3.3  Let  Pn{x)  be  the  n-th  normalized  Legendre  polynomial  (defined  in  (IS)). , 
Then  for  any  real  a, 

j^^e^‘^^K{x)dx 

QQ  rX  ___  ^  rl 

=  X^'' Pn{x)  dx i  '^fik  Pn{x)  dx.  (39) 

k=ko  k=ko 


where 


2k 

(j2k+l 

^  (2A:  +  1)!’ 

kQ  =  [n/2\ . 


(40) 

(41) 

(42) 


Furthermore,  for  all  integer  m>  [e-  |a|J  +  1, 


(43) 
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In  particular,  if 

^  >2  (Le-  |a|J  +1)  , 

then 


(44) 


(45) 


Proof.  The  formula  (39)  follows  immediately  from  Lemma  2.2  and  Taylor’s  expansion 
of  In  order  to  prove  (43),  we  assume  that  m  is  an  integer  such  that 


m  >  [e  •  |a|J  +  1 . 

Introducing  the  notation 

Rm=  Y^Olk  P„(rc)  dx  +  iY^  Pk  /  Pn{x)  dx, 

k=m  “'-1  •'-1 


(46) 


(47) 


k=m 


we  immediately  observe  that,  due  to  Lemma  2.2  and  the  triangle  inequality, 

lai'^  rr 


|i2m|  < 

00 

E 

fc=2m 

< 

00 

E 

fc=2m 

Since  (46)  implies  that 

k\  y  k  +  1 

Ik 


a 


1^ 

k\ 


a\  1  1 

<  TT-  <  77-  < 


(48) 


2m  4-  k  2m  2e  2  ’ 
for  all  integer  m,  fc  >  0,  we  rewrite  (48)  as 


(49) 


I-Rmi  < 


(2m)! 

I^|2m 


<  2 


(50) 


(2m) !  ’ 

and  obtain  (43)  immediately  using  Stirling’s  formula.  Finally,  we  obtain  (45)  by  choosing 
m  =  (e  •  |a|J  +  1 .  (51) 

□ 


11 


Theorem  3.4  Let  be  the  m-th  prolate  spheroidal  function  with  band  limit  c,  let 

Pk{x)  be  the  k-th  normalized  Legendre  polynomial  (defined  in  (13)),  and  let  Xm  be  the 
eigenvalue  which  correspands  to  iprni^)  (qs  in  Theorem  2.4)-  Then  for  all  integer  m  >  0 
and  all  real  positive  c,  if 

A:>2  ([e-cJ  +  1)  ,  (52) 

then 

j^^'iprn{x)P'k{x)  dx  <  ~  .  (53) 

Moreover,  given  any  £>  0,  if 

/c  >  2  ( [e  •  cj  +  1)  +  logs  +  log2  (3^) .  (54) 

then 

y  ^'il)m{x)Ti,{x)  dx  <  £  .  (55) 

Proof.  Obviously 

^j^^'(l;m{x)T'k{x)  dx 

^  /  lV'm(2:)|*|/  e^^^P'k{t)dtdx.  (56) 

Introducing  the  notation 

a  =  cx,  (57) 


and  remembering  that 

J  ^\i)rn{x)\dx  =  l,  (58) 

we  observe  that  the  combination  of  (56),  (57),  (58),  and  Lemma  3.3  implies  that 

^i’mix)  Pk{x)  dx^ 

<  io'(D 


1  n 


|A^I  V2;  • 

Substituting  (54)  into  (53),  we  immediately  see  (55). 
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4  Numerical  Evaluation  of  Prolate  Spheroidal  Wave- 
functions 

Both  the  classical  Bouwkamp  algorithm  (see,  for  example,  [1])  for  the  evaluation  of  the 
functions  ipj,  and  the  algorithm  presented  in  this  paper  for  the  same  task,  are  based  on 
the  expression  of  those  functions  as  a  Legendre  series  of  the  form 


=  E  <^kPk(x)-, 


(60) 


k=(i 


since  the  functions  ij^j  are  smooth,  the  coefficients  ak  decay  superalgebraically  (with 
bounds  for  that  decay  being  given  in  Theorem  3.4).  Substituting  (60)  into  (24),  and 
using  (8)  and  (11),  we  obtain  the  well-known  three-term  recursion 


■  ak+2  + 


(A:  -h  2)(A;  -t- 1) 

(2k  -I-  3)(2k  -(-  5) 

+  1)  +  -In  • 


k(k-l) 


{2k  +  3)(2k-l) 

•  (P  •  ak-2  =  0. 


c^k  + 


(61) 


(2k  -  3)(2A:  -  1) 

Combining  (61)  with  (13),  we  obtain  the  three-term  recursion 

(k  -f  2)(A:  4- 1) 


(2k  4-  3)^(2k  +  5)(2k  +  l) 

A.,  2A:(fc4-l)-l  2 

fe(^  +  l)-k  vwo,-  nx-g  -Xj 


(2k  +  3)(2k-l) 
k(k  -  1)  2 


Pi  + 


(62) 


(2k  -  l)y/(2k-3)(2k+l) 

for  the  coefficients  of  the  expansion 

00 

=  E  Pk-^k(^)\ 


c"-/5L2  =  0 


(63) 


k=0 


for  each  j  =  0, 1, 2, ... ,  we  will  denote  by  the  vector  in  P  defined  by  the  formula 

p^  =  (pipi,pi...).  (64) 


The  following  theorem  restates  the  recursion  (62)  in  a  slightly  different  form. 


Theorem  4.1  The  coefficients  Xi  are  the  eigenvalues  and  the  vectors  are  the  corre¬ 
sponding  eigenvectors  of  the  operator  represented  by  the  symmetric  matrix  A 

given  by  the  formulae  ‘ 


,  ,,,  2A:(A:  + 1)  - 1  , 

A,,,  -  M<:+l)+(2t+3)(24_l)-"’ 

(65) 

(A:  +  2)(A: -i- 1)  2 

^k,h+2  —  ! -  *  c  , 

{2k4-^)^J{2k-\-l){2k-\-5) 

(66) 

(k  +  2){k-\-l)  2 

■A:+2,fc  =  / -  •  C  , 

{2k  +  3)^{2k-i-l){2k  +  5) 

(67) 

for  all  k  =  0,1,2,. . with  the  remainder  of  the  entries  of  the  matrix  being  zero. 

In  other  words,  the  recursion  (62)  can  be  rewritten  in  the  form 

iA-xj^I)m^0,  (68) 

where  A  is  separable  into  two  symmetric '  tridiagonal  matrices  and  A^dd,  the  first 
consisting  of  the  elements  of  A  with  even-numbered  rows  and  columns  and  the  second 
consisting  of  the  elements  of  A  with  odd-numbered  rows  and  columns.  While  these  two 
matrices  are  infinite,  and  their  entries  do  not  decay  much  with  increasing  row  or  column 
number,  the  eigenvectors  of  interest  (those  corresponding  to  the  first  m  prolate 
spheroidal  functions)  lie  almost  entirely  in  the  leading  rows  and  columns  of  the  matrices 
(as  shown  by  Theorem  3.4).  Thus  the  evaluation  of  prolate  spheroidal  functions  can  be 
performed  by  the  following  procedure: 

•  1.  Generate  the  leading  k  rows  and  columns  of  A,  where  k  is  given  by  (54). 

•  2.  Split  the  generated  portion  of  A  into  and  A^dd,  and  use  a  solver  for  the 
symmetric  tridiagonal  eigenproblem  (such  as  that  in  LAPACK)  to  compute  their 
eigenvectors  {/3^  and  eigenvalues  {Xj}- 

•  3.  Use  the  obtained  values  of  the  coefficients  Pi,  P{,  . . .  in  the  expansion  (63)  to 

evaluate  the  function  ■j/’j  at  arbitrary  points  on  the  interval  [—1,1]- 

Obviously  steps  1  and  2  can  be  performed  as  a  precomputation,  for  any  given  value  of 
c.  As  a  numerical  diagonalization  of  a  positive  definite  tridiagonal  matrix  with  well- 
separated  eigenvalues,  this  precomputation  stage  is  numerically  robust  and  efficient, 
requiring  0(cm)  operations  to  construct  the  Legendre  expansions  of  the  form  (64)  for  the 
first  m  prolate  spheroidal  functions;  each  subsequent  evaluation  of  a  prolate  spheroidal 
function  takes  0(c)  operations. 
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4.1  Numerical  Evaluation  of  Eigenvalues 

Although  the  above  algorithm  for  the  evaluation  of  prolate  spheroidal  wave  functions  also 
produces  the  eigenvalues  {xj}  of  the  differential  operator  (24),  it  does  not  produce  the 
eigenvalues  {Aj}  of  the  integral  operator  Fc  (defined  in  (17)).  Some  of  those  eigenvalues 
can  be  computed  using  the  formula 

Xjipj{x)  =  j  dt,  (69) 

evaluating  the  integral  on  the  right  hand  side  numerically;  however,  that  evaluation 
obviously  has  a  condition  number  of  about  l/\j,  and  is  thus  inappropriate  for  computing 
small  Xj.  A  well-conditioned  procedure  is  as  follows: 

•  1.  Use  (69)  to  calculate  Aq,  evaluating  the  right  hand  side  numerically,  and  with 
rr  =  0  (so  that  'ipo(x)  is  not  small). 

•  2.  Use  the  calculated  Aq,  together  with  Corollary  3.2,  to  compute  the  absolute  val¬ 
ues  [Ajl,  for  j  =  1,2,  computing  each  |Aj|  from  |Aj_i|  (and  again,  evaluating 
the  required  integrals  numerically). 

•  3.  Use  the  fact  that  Xj  =  i-^lAj]  (see  Theorem  2.4)  to  finish  the  computation. 

5  Quadratures  for  Band-Limited  Functions 

Since  the  prolate  spheroidal  wave  functions  •  •  •  >  V’nj  •  •  •  constitute  a  complete  or¬ 

thonormal  basis  in  L^[— 1,1]  (see  Theorem  2.4), 

=  E  (/' M*),  (70) 

j=Q  V-1  / 

for  all  x,t  £  [—1, 1];  substituting  (18)  into  (70)  yields 

C» 

gicxt  ^  ^ 

Thus  if  a  quadrature  integrates  exactly  the  first  n  eigenfunctions,  that  is,  if 


for  all  j  =  0, 1, . . .  ,n— 1,  then  the  error  of  the  quadrature  when  applied  to  a  function 
f{x)  =  with  a  G  [—1, 1],  is  given  by 

m  ,  -i 

T  Wke^‘^^>‘  -  /  dx 

m  I  OQ  \  rl  #  \ 

=  Z!  I]  V’j (<i)  i’jM]-  Z  (^)  V’j (2;) 

k=l  \j=Q  I  •'-1  \j=Q  J 


m/00  \  .1  00  \ 

=  Z  (  Z  (^)  %  (^Jk)  j  -  J_^  (  Z  (“)  V’i  (^)  j  dx.  (73) 

Due  to  the  orthonormality  of  the  functions  {ipj}, 

00 

E  (74) 


j=zn 


Wr. 


Prom  (74),  it  is  obvious  that  the  error  of  integration  (73)  is  of  roughly  the  same  mag¬ 
nitude  as  An,  provided  that  n  is  in  the  range  where  the  eigenvalues  {Aj}  are  decreasing 
exponentially  (as  is  the  case  for  quadratures  of  any  useful  accuracy;  see  Theorem  2.5) 
and  provided  in  addition  that  the  weights  {ly*}  are  not  large. 

Now,  the  existence  of  an  n/2-point  quadrature  that  is  exact  for  the  first  n  Prolate 
Spheriodal  Wave  functions  follows  from  the  combination  of  Theorems  2.1,  2.4;  an  al¬ 
gorithm  for  the  numerical  evaluation  of  nodes  and  weights  of  such  quadratures  can  be 
found  in  [2].  An  alternative  procedure  for  the  construction  of  quadrature  formulae  for 
band-limited  functions  (leading  to  slightly  different  nodes  and  weights)  is  described  in 
the  following  section;  a  numerical  comparison  of  the  two  can  be  found  in  Section  8  below. 


Remark  5.1  The  above  text  considers  only  the  error  of  integration  of  a  single  exponen¬ 
tial.  For  a  band-limited  function  g  :  [— 1, 1]  — >  C  given  by  the  formula 

9ix)=  (75) 

for  some  function  G  :  [—1, 1]  C,  the  error  is  obviously  bounded  by  the  formula 


Y^Wkg{xk)-l  g{x)  dx 

fc=i 


(76) 


where  e  is  the  maximum  error  of  integration  (73)  of  a  single  exponential,  for  any  t  G 
[—1, 1].  While  ||G1|  might  be  much  larger  than  (as  it  is  if,  for  instance,  g  =  V’so-n), 

if  the  same  equation  (75)  is  used  to  extend  g  to  the  rest  of  the  real  line,  then  by  Parseval’s 
formula  ||G1|  =  l|5l|(-oo,oo);  that  is  to  say,  although  the  error  of  such  a  quadrature  when 
applied  to  a  band-limited  function  is  not  bounded  proportional  to  the  norm  of  that 
function  on  the  interval  of  integration,  it  is  bounded  proportional  to  the  norm  of  that 
function  on  the  entire  real  line. 


6  Quadrature  Nodes  from  Roots  of  Prolate  Func¬ 
tions 

An  alternative  to  the  approach  of  the  previous  section  is  to  use  roots  of  appropriate 
prolate  spheroidal  wave  functions  as  quadrature  nodes,  with  the  weights  determined  via 
the  procedure  described  in  Remark  2.1.  The  following  theorems  provide  a  basis  for  this; 
numerically  (see  Section  8)  the  resulting  quadrature  nodes  tend  to  be  inferior  to  those 
produced  by  the  optimization  scheme  of  [13,  25,  2];  however,  they  are  useful  as  starting 
points  for  that  scheme,  or  as  somewhat  less  efficient  nodes  which  can  be  computed  much 
more  quickly. 

6.1  Euclid  Division  Algorithm  for  Band-Limited  Functions 

The  following  two  theorems  constitute  a  straightforward  extension  to  band-limited  func¬ 
tions  of  Euclid’s  division  algorithm  for  polynomials’.  Their  proofs  are  quite  simple,  and 
are  provided  here  for  completeness,  since  the  author  failed  to  find  them  in  the  literature. 

Theorem  6.1  Suppose  that  cr,  :  [0, 1]  C  are  a  pair  of  c^— functions  such  that 


y'(  1)7^0,  (77) 

c  is  a  positive  real  number,  and  the  functions  /,  p  are  defined  by  the  formulae 

f{x)=  f' (78) 
Jo 

p{x)  =  f  (p{t)  6*“*  dt.  (79) 

Jo 

Then  there  exist  two  c^ -functions  :  [0, 1]  — >■  (D  such  that 

f{x)  =  p{x)  q(x)  +  r{x)  (80) 

for  all  a;  €  R,  with  the  functions  q,  r  :  [0, 1]  — >  R  defined  by  the  formulae 

q(x)  =  [  T]{t)  dt,  (81) 

Jo 

r(x)  =  /'e(i)e'“‘ dt.  (82) 

Jo 
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Figure  1:  The  split  of  integration  range  that  yields  (85) 


Proof. 

Obviously,  for  any  functions  p,q  given  by  (79),  (81), 

p{x)  q{x)  =  f  (p(t)  dt  •  f  7]{t)  dr 
Jo  Jo 


—  f  f  ¥’(t)  ??(v)  drdt. 

Jo  Jo 

Defining  the  new  independent  variable  u  by  the  formula 

U  =  t  +  T, 

we  rewrite  (83)  as 

p{x)q(x)  =  [  [  (p{u—T)ri{T)  drdu 

Jo  Jo 

+  [  [  (p{u—t)p(t)  drdu 

J\  Ju  —  \ 

(see  Figure  1).  Substituting  (78),  (82),  and  (85)  into  (80),  we  get 

[  f  (p(u—t)p(t)  drdu 
Jo  Jo 

4-  f  [  ip{u—T)p{T)  dr  du-^  f  ^(t)e*“* 
Ji  Ju-i  Jo 

/•1/2  /■! 

=  /  alt)  dt  +  /  a{t)  dt. 

Jo  Jl/2 


dt 


(83) 

(84) 


(85) 


(86) 


Due  to  the  well  known  uniqueness  of  the  Fourier  Transform,  (86)  is  equivalent  to  two 
independent  equations; 


[  f  (p(u—r)p{T)  dr  du+  f  ^{t)e^^^dt 

Jo  Jo  Jo 

ri/2 

=  /  a(t)  dt, 

Jo 


(87) 
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f  ^tcux  f  ip(u—r)T]{r)drdu=  f  a{t)  dt.  (88) 

Jl  Ju—l  Jl/2 

Now,  we  observe  that  (88)  does  not  contain  and  use  it  to  obtain  an  expression  for  rj 
as  a  function  of  cp,  a.  After  that,  we  will  view  (87)  as  an  expression  for  ^  via  (p,  a,  rj. 
From  (88)  and  the  uniqueness  of  the  Fourier  Transform,  we  obtain 

[  (p{u-t)  r){r)  dr  =  a(^),  (89) 

Ju—l  L 

for  all  tt  6  [1, 2].  Introducing  the  new  variable  v  via  the  formula 

V  =  u-l,  (90) 


we  convert  (89)  into 


L 


1  w  +  1 

ip{v+l-T)t)(T)  dr  =  <r(— — ), 


(91) 


which  is  a  Volterra  equation  of  the  first  kind  with  respect  to  differentiating  (91)  with 
respect  to  u,  we  get 

-¥’(1)  Viv)  +  ^  (p'iv+l-r)  7]{t)  dT  =  ^  (92) 

which  is  a  Volterra  equation  of  the  second  kind.  Now,  the  existence  and  uniqueness 
of  the  solution  of  (92)  (and,  therefore,  of  (89)  and  (88))  follows  from  Theorem  2.3  of 
Section  2. 

With  Tj  defined  as  the  solution  of  (89),  we  use  (87)  together  with  the  uniqueness  of 
the  Fourier  Transform,  to  finally  obtain 

^(w)  =  ^  v{u-r)r]{T)  dr,  (93) 


for  all  u  €  [0, 1]. 

The  following  theorem  is  a  consequence  of  the  preceding  one. 


□ 


Theorem  6.2  Suppose  that  a,(p  :  [—1,1]  — >•  C  are  a  pair  of  c^— functions  such  that 
<p{—l)  ^  0,  (;^j(l)  7^  0,  c  is  a  positive  real  number,  and  the  functions  f,p  are  defined  by 
the  formulae 


f{x)  =  J  ^  a{t)  dt, 
p{x)  =  J  ip{t)  e'“^^  dt. 


(94) 

(95) 
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Then  there  exist  two  -functions  r],^  :  [—1,1]  C  such  that 
f(x)  =  p{xlq{x)  +  r{x) 

for  all  a;  €  R,  with  the  functions  q,  r  :  [—1, 1]  — >■  R  defined  by  the  formulae 
q(x)  =  dt, 

r(l)  =  r  dt. 

Proof. 

Defining  the  functions  /+,  f-,p+,  P-,  by  the  formulae 

/+(x)  =  l\(t)^^‘dt, 

Jo 

/-W=  dt, 

p+{x)  =  f  -Pit)  dt, 

Jo 

p-ix)  -  j  ^(p{t)e'‘^^  dt, 

we  observe  that  for  all  x  €  R\ 

fix)  =  U{x)+f-ix), 
p(x)  =p+(x}  -hp-(x). 

Due  to  Theorem  6.1,  there  exist  such  r]+,  rj-,  ^+,  that 
/+  (x)  =  p+  (x)  q+  (x)  +  r+  (x) , 

/_(x)  =  p_(x)  9_(x)  +  r_(x), 
with  the  functions  q^,q-,r+,r-  defined  by  the  formulae 

q+(x)  =  dt, 

Jo 

s_(i)  =  y  )?-(«)  e*"'  dt, 


(96) 

(97) 

(98) 

(99) 
(100) 
(101) 
(102) 

(103) 

(104) 

(105) 

(106) 

(107) 

(108) 
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^+(2;)  =  /  ?+(i)e*“‘  dt, 

Jo 

(109) 

r.{x)  = 

(110) 

Now,  defining  q,  by  the  formula 

q{x)  =q4x)  +  q+{x) 

(111) 

for  all  X  G  [—1, 1],  we  have 

pix)q{x)  =  {p.(x)  +  p+(x))  ■  (q-ix)  +  q+{x)) 

=  p+{x)q+{x)+p_{x)q-{x)+p-{x)q+(x)+p+{x)q.{x),  (112) 

and  we  define  r{x)  by  the  obvious  formula 

r{x)  =  r-{x)  +  r+(x)  -  (p-{x)  q+{x)  +p+(x)  q-{x)).  (113) 

□ 


6.2  Quadrature  nodes  from  the  division  theorem 

In  much  the  same  way  that  the  division  theorem  for  polynomials  can  be  used  to  provide 
a  constructive  proof  of  Gaussian  quadratures,  Theorem  6.2  provides  a  method  of  con¬ 
structing  generalized  Gaussian  quadratures  for  band-limited  functions.  The  method  is 
as  follows. 

To  construct  a  quadrature  for  functions  of  a  bandwidth  2c,  prolate  spheroidal  wave 
functions  corresponding  to  bandwidth  c  are  used.  (Thus  the  eigenvalues  {Aj}  and  eigen¬ 
functions  {i/jj}  are  in  this  section,  as  elsewhere  in  the  paper,  those  corresponding  to 
bandwidth  c).  The  following  theorem  provides  a  bound  of  the  error  of  a  quadrature 
whose  nodes  are  the  roots  of  the  n’th  prolate  function  when  applied  to  a  function 
/  which  satisfies  the  conditions  of  the  division  theorem,  in  terms  of  the  norms  of  the 
quotient  and  remainder  of  /  divided  by  V’n: 

Theorem  6.3  Suppose  that  xi,X2, . . .  ,Xn  €  R  are  the  roots  of  ipn-  Let  the  numbers 
wi,W2, . . . ,  rcn  G  R.  be  such  that 

^  rl 

Y^Wk'4>j{xk)  =  j  ipj{x)  dx,  (114) 

fc=i 
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for  all  j  =  0, 1, . . . ,  n  —  1.  Then  for  any  function  f  :  [^1, 1]  — >•  C  which  satisfies  the 
conditions  of  Theorem  6.2, 

n  ■'  1 

Wkf{xk)  -  /  fix)  dx 

k=l 

00  /  m  \ 

<  iA„i  ■  ii»ii  +  lien  ■  r  iAii  ■  WiiiL  •  2 + x:  iKii  .  (115) 

j=n  \  k=l  ) 

where  the  functions  :  [—1, 1]  C  are  as  defined  in  Theorem  6.2. 

Proof.  Since  /  satisfies  the  conditions  of  Theorem  6.2,  there  exist  functions  g,r  : 
[-1,1]  1,  defined  by  (97), (98)  such  that 

fix)  =  tfjnix)  qix)  +  rix).  (116) 


Then,  defining  the  error  of  integration  Ef  for  the  function  /  by 


Ef  = 


^  ^  pi 

Y  ^kfiXk)  -  /  fix)  dx 

k=l 


(117) 


we  have 


Ef  = 


< 


Y  'Wk  i'ipnixk)  qixk)  +  r(a;*))  -  /  i-ipnix)  qix)  +  r(a;))  dx 

k=l 

n  -1 

Y  Wk  IpniXk)  qiXk)  -  /  ipnix)  qix)  dx 

”  rl 

Y  '^krixk)  -  /  rix)  dx 
1 _ 1  1 


Since  the  nodes  are  the  roots  of  if., 


715 


Y'^ki^niXk)qiXk)  =  0. 

k=l 


Thus 


Now 


<  I  /  ipnix)qix)  dx  +  Y'^k'rixk)  -  [  rix)  dx 
U-i  J-i 

j_^if nix)  qix)dx  =  J  V’n(a^)  f  v{^) 

=  j  T](t)  J  Tpn(x)  6*“®*  dx  dt 

=  j  n(t) Ki’A*) 


(118) 


(119) 


(120) 


(121) 
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Using  the  Cauchy-Schwartz  inequality  and  the  fact  that  the  function  'ipn  has  unit  norm, 
we  get  from  (121)  that 


Also, 


q{x)  dx 


(122) 


Wkr{xk)  -  r{x)  dx 

=  E  (/_' «‘) *)  -  (£  f  (*)  *) 

=  dxjdt.  (123) 

Substituting  (73)  into  (123),  and  using  the  Cauchy-Schwartz  inequality,  we  get 

n  -1 

Y^krixk)-  r(x)  dx 

k=i 

-1  /  m  /  oo 

•'-1  \k=l  \j=n 

OO  /  \ 

<  2  +  EIKII  ■  (124) 

j—n  \  A:=l  / 

Combining  (120),  (122),  and  (124),  we  get 

OO  /  m  \ 

E,  <  |A„|-I|.)||  +  IKII-El^ll-lll^’illl-  2  +  ElKII  ■  (125) 

j=n  \  fc=l  / 

□ 


Remark  6.1  The  use  of  Theorem  6.3  for  the  construction  of  quadrature  rules  for  band- 
limited  functions  depends  on  the  fact  that  the  norms  of  the  band-limited  functions  q 
and  r  in  (116)  are  not  large,  compared  to  the  norm  of  /  (both  sets  of  norms  being  on 
[—00,00]).  Such  estimates  have  been  obtained  for  all  n  >  2c/7r  -1-  lOlog(c).  The  proofs 
are  quite  involved,  and  will  be  reported  at  a  later  date.  In  this  paper,  we  demonstrate 
the  performance  of  the  obtained  quadrature  formulae  numerically  (see  Section  8  below). 


Remark  6.2  It  is  natural  to  view  (116)  as  an  analogue  for  band-limited  functions  of 
the  Euclid  division  theorem  for  polynomials.  However,  there  are  certain  differences.  In 
particular,  Theorem  6.1  admits  extensions  to  band-limited  functions  of  several  variables, 
while  the  classical  Euclid  algorithm  does  not.  Such  extensions  (together  with  several 
applications)  will  be  reported  at  a  later  date. 

7  Interpolation  via  Prolate  Spheroidal  Wavefunctions 

Interpolation  is  usually  performed  by  the  following  general  procedure;  assuming  that  the 
function  /  :  [a,  6]  C  to  be  interpolated  is  given  by  the  formula 

f{x)  =  Ci0i(2:)  -h  C2(t)2ix)  +  .  .  .  +  Cn(l>nix),  (126) 

where  <l)i,<f>2,  ■  ■  •  ,  <f>n  ■  [<3^5^]  are  a  fixed  sequence  of  functions  (often  polynomials), 
solve  an  n  X  n  linear  system  to  determine  the  coefficients  Ci,  C2, . . . ,  c„  from  the  values  of 
/  at  the  n  interpolation  nodes,  then  use  (126)  to  evaluate  /  wherever  needed.  As  is  well 
known,  if  /  is  well-approximated  by  a  linear  combination  of  the  interpolation  functions, 
and  if  the  linear  system  to  be  solved  is  well-conditioned,  then  this  procedure  is  accurate. 

As  shown  in  Section  5  in  the  context  of  quadratures,  a  linear  combination  of  the  first 
n  prolate  spheroidal  functions  ^O) V’l) •  •  •  for  a  band  limit  c  can  provide  a  good 
approximation  to  functions  of  the  form  with  f  €  [-1, 1]  (see  (71,74));  in  the  regime 
where  the  accuracy  is  numerically  useful,  the  error  is  of  the  same  order  of  magnitude  as 
|A„|.  This,  in  turn,  shows  that  they  provide  a  good  approximation  (in  the  same  sense  as 
in  Remark  5.1)  to  any  band-limited  function  of  band  limit  c.  Thus,  if  ■00)  V'l  ?  ■  •  •  >  V’n-i  are 
used  as  the  interpolation  functions  in  this  procedure,  they  can  be  expected  to  yield  an 
accurate  interpolation  scheme  for  band-limited  functions,  provided  that  the  matrix  to  be 
inverted  is  well-conditioned.  The  following  theorem  shows  that  if  the  interpolation  nodes 
are  chosen  to  be  quadrature  nodes  accurate  up  to  twice  the  bandwidth  of  interpolation, 
with  the  quadrature  formula  being  accurate  to  more  than  twice  as  many  digits  as  the 
interpolation  formula  is  to  be  accurate  to,  then  the  matrix  inverted  in  the  procedure  is 
close  to  being  a  scaled  version  of  an  orthogonal  matrix. 

Theorem  7.1  Suppose  the  numbers  W\,u)2,  G  R  and  xi,X2, .  ■  ■  ,Xn  ^  ^  are  such 

that 


for  all  a  €  [—1, 1],  and  for  some  c  >  0.  Let  the  matrix  A  be  given  by  the  formula 

[  A{^l)  ■■■  \ 

^’0(2:2)  ^1(2:2)  •••  ■^71-1(3:2) 


(128) 


\'tpo{Xn)  i>\{Xn)  •••  i’n-lM  J 

let  the  matrix  W  be  the  diagonal  matrix  whose  diagonal  entries  are  Wi,W2, . .  .  ,Wn,  and 
let  the  matrix  E  =  [ejk]  be  given  by  the  formula 

E  =  I-  A*WA. 

Then 

2e 


\^ik\  < 

Proof.  Clearly 


Aj_iA 


fe-i 


(129) 

(130) 


(131) 


Cjk  =  Sjk  i’j-lixi)  -ifk-iixi), 

i=i 

where  Sij  is  the  Kronecker  delta  function.  Using  (18),  this  becomes 

ejk  =  Sjk-J2'^i-  (t~  [  e~^'^^^ifj-i{t)dt 

i=i  ^ 

\Ak-i  J-1  j 

=  ^jk-= — \ - f  f  'il^j-i{t)'ipk-i{T)Y^wie~^^‘^e"‘^^^dtdT.  (132) 

Aj—iAk—l  i—i 

Using  (127),  this  becomes 

ejk  —  djk~=  r  f  f  'tpj-i{t)  ipk-i{T)  (133) 

■'-1  •'-i  /  •  \ 

•  ij  e  —  fe{t+r)]  dtdr, 

where  /e :  [— 2, 2]  — >  C  is  a  function  which  satisfies  the  relation 

\feix)\<e,  (134) 

for  all  X  €  [—2, 2].  Thus 

ejk  =  ^jk--F — \ - f  f  i)j-i{t)ifk-\{r)  !  ds dt dr 

Aj—\Ak—\  J J J—l 

+= — \ — [  f  ipj-i{t)xpk-i{r)  fe(t  +  T)  dtdr  (135) 

•'-1  •'-1 
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Using  (18),  this  becomes 

ejk  =  Sjk- J^^  'ipj-i{s)ipk-i{s)  ds 

+Y~\ —  /  '>Pk-i{'r)  f  ijjj-iit)  feit  +  t)  dtdr. 

Aj-lAk-l  •> -1-  J—l 

Due  to  the  orthonormality  of  the  functions  {ipj},  this  becomes 

1  /•!  1-1 

ejk  =  Y — r -  /  ^fc-i(r)  /  fe{t  +  r)  dtdr. 

Aj—iAfe-i  J-1  J—1 

Using  the  Cauchy-Schwartz  inequality,  this  becomes 

1 


\ejk\  ^ 


< 


1 


^j-l^k-l 

1 


llV-fc-ill  ]j fe{t  +  r)  dt 


dr 


2e 


A,_iAfc_i 


+  dt  dr 


(136) 


(137) 


(138) 


□ 


From  inspection  of  Theorem  2.5,  it  can  easily  be  seen  that  the  number  N  of  eigenval¬ 
ues  needed  for  a  bandwidth  of  2c  and  an  accuracy  of  is  roughly  twice  the  number  of 
eigenvalues  needed  for  a  bandwidth  of  c  and  an  accuracy  of  e.  Thus  a  generalized  Gaus¬ 
sian  quadrature  for  a  bandwidth  2c  and  an  accuracy  has  roughly  the  same  number 
of  nodes  as  are  needed  for  interpolation  of  accuracy  e.  In  our  numerical  experiments, 
this  correspondence  was  found  to  be  much  closer  than  the  rough  bounds  in  Theorem  2.5 
indicate;  in  the  results  tabulated  in  Section  8,  the  number  of  nodes  for  an  interpolation 
formula  of  a  desired  accuracy  e  was  always  chosen  to  be  the  number  of  quadrature  nodes 
for  a  desired  accuracy  for  twice  the  band  limit  (that  number,  in  turn,  being  chosen 
as  indicated  in  Section  5);  the  correspondence  between  the  desired  accuracy  and  the 
experimentally  measured  maximum  error  can  be  seen  in  Tables  3  and  4. 

The  coefficients  ci,C2,  ...,c„  produced  by  this  interpolation  procedure  (see  (126)) 
can,  of  course,  just  as  easily  be  used  for  evaluating  derivatives  or  indefinite  integrals  of 
the  interpolated  function,  as  they  can  for  computing  the  function  itself. 
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8  Numerical  Results 


The  algorithms  of  Sections  5-7  have  been  implemented  in  double  precision  (64-bit  floating 
point)  arithmetic,  with  results  shown  in  Tables  1-4.  Tables  1  and  2  show  the  perfor¬ 
mance  of  quadrature  nodes  produced  by  the  schemes  of  Sections  5  and  6,  when  used  as 
quadrature  nodes;  Tables  3  and  4  show  their  performance  when  used  as  interpolation 
nodes.  These  are  not  actually  the  same  sets  of  nodes;  even  with  the  bandwidth  c  for  in¬ 
terpolation  being  half  of  the  bandwidth  for  quadrature  (as  it  is  in  the  tables),  more  nodes 
are  needed  to  achieve  a  given  accuracy  of  interpolation  than  are  needed  to  achieve  a  given 
accuracy  of  quadrature,  as  can  be  seen  by  comparing  the  number  of  nodes  (printed  in 
the  column  labeled  n  in  each  table).  The  error  figures  in  the  tables  are  approximations 
of  the  maximum  error  of  interpolation  or  of  the  quadrature,  when  applied  to  functions 
of  the  form  cos(ax)  and  sin(arc),  with  0  <  a  <  c;  they  were  computed  by  measuring  the 
error  at  a  large  number  of  points  in  a  (for  interpolation,  in  both  a  and  x),  including  the 
extremes.  The  column  labeled  “Roots”  contains  the  errors  for  the  nodes  produced  by 
the  scheme  of  Section  6;  the  column  labeled  “Refined”  contains  the  errors  after  those 
nodes,  used  as  a  starting  point,  have  been  run  through  the  scheme  of  Section  5.  The 
variable  e  which  appears  in  the  tables  is  the  requested  accuracy,  used  to  determine  the 
number  of  nodes  in  the  ways  described  in  Sections  5  and  7. 

Also  tabulated  are  the  numbers  of  Legendre  nodes  required  to  achieve  the  same 
accuracy  s  using  polynomial  interpolation  or  quadrature  schemes.  Since  Chebyshev 
nodes  are  generally  known  to  be  superior  for  interpolation,  for  that  case  the  numbers  of 
Chebyshev  nodes  required  to  achieve  the  same  accuracy  are  also  tabulated. 

Figure  2  contains  the  maximum  norm  of  the  derivative  of  each  prolate  function  i}j{x), 
for  c  =  200  and  x  G  [-1,1],  as  a  function  of  also  graphed,  for  comparison,  is  the 
maximum  norm  of  the  derivative  of  each  normalized  Legendre  polynomial  Pj{x)  over 
the  same  range;  and  graphed  below,  on  the  same  horizontal  scale,  are  the  norms  of  the 
eigenvalues  \j.  The  graph  shows  that,  for  this  value  of  c,  computing  the  derivatives  of 
a  function  given  by  a  prolate  series  is  a  better-conditioned  operation  than  computing 
the  derivatives  of  a  function  given  by  a  Legendre  series  of  the  same  number  of  terms. 
(Obviously,  if  the  number  of  terms  can  also  be  reduced,  as  in  the  situations  of  Tables  1- 
4,  there  is  a  further  improvement  in  the  condition  number.)  The  same  general  pattern 
of  behavior  is  exhibited  for  other  values  of  c;  as  c  approaches  zero  (and  the  prolate 
functions  approach  the  Legendre  polynomials),  the  value  of  j  at  which  the  maximum 
norm  of  the  derivative  rises  sharply  also  approaches  zero  (as  is  to  be  expected,  since  for 
c  =  0  the  prolate  functions  reduce  to  Legendre  polynomials).  Finally,  Tables  5  and  6 
contain  samples  of  quadrature  weights  and  nodes. 

Remark  8.1  In  this  paper,  detailed  discussion  of  issues  encountered  in  the  implemen¬ 
tation  of  numerical  algorithms  has  been  deliberately  avoided,  as  well  as  any  discussion  of 


CPU  time  requirements,  memory  requirements,  etc.  Thus,  we  limit  ourselves  to  observ¬ 
ing  that  all  algorithms  have  been  implemented  in  FORTRAN,  that  with  the  exception  of 
the  procedure  for  the  evaluation  of  Prolate  Spheroidal  Wave  functions  described  in  Sec¬ 
tion  4,  we  have  not  designed  or  implemented  any  new  or  original  numerical  algorithms, 
and  that  the  procedure  of  Section  4  consists  of  applying  standard  tools  of  numerical 
analysis  (diagonalization  of  a  tridiagonal  matrix)  to  the  well-known  recursion  (61).  The 
resulting  algorithm  for  the  evaluation  of  prolate  spheroidal  wave  functions  has  the  CPU 
time  requirements  proportional  to  (?,  with  a  fairly  large  proportionality  constant.  The 
procedure  of  [2],  when  applied  to  the  system  of  functions  ^o>  •  •  •  j  ^2n+i  requires  order 

T?  operations,  also  with  a  fairly  large  proportionality  constant.  On  the  other  hand,  the 
cost  of  finding  all  roots  n  of  the  function  ipn  lying  on  the  interval  [—1, 1]  is  proportional 
to  n,  and  the  proportionality  constant  is  not  large.  The  largest  c  we  have  dealt  with  in 
our  experiments  was  about  6000,  with  resulting  quadratures  having  about  1900  nodes. 
In  this  regime,  the  construction  of  the  quadrature  (both  nodes  and  weights)  took  several 
minutes  on  the  300-m^g^op  SUN  workstation;  while  there  are  fairly  obvious  ways  to 
reduce  the  cost  of  the  calculation  (both  in  terms  of  asymptotic  CPU  time  requirements 
and  in  terms  of  associated  proportionality  constants)  we  have  made  no  effort  to  do  so. 

The  following  observations  can  be  made  from  the  examples  presented  in  this  section, 
and  from  the  more  extensive  tests  performed  by  the  authors. 

1.  When  the  nodes  obtained  via  the  algorithm  of  [2]  are  used  for  the  integration  of  band- 
limited  functions,  the  resulting  quadrature  rules  are  significantly  more  accurate  than  the 
quadratures  obtained  from  the  nodes  of  appropriately  chosen  prolate  functions;  however, 
the  difference  between  the  numbers  of  nodes  required  by  the  two  approaches  to  obtain 
a  prescribed  precision  is  not  large.  When  the  nodes  obtained  via  the  two  approaches  are 
used  for  the  interpolation  (as  opposed  to  the  integration)  of  band-limited  functions,  the 
performances  of  the  two  are  virtually  identical. 

2.  For  large  c,  the  number  of  nodes  required  by  a  quadrature  rule  for  the  integration 
of  band-limited  functions  with  the  band-limit  c  is  close  to  the  dependence  on  the 
required  precision  of  integration  is  weak  (as  one  would  expect  from  Theorem  2.5  and 
subsequent  developments). 

3.  The  numbers  of  nodes  required  by  our  quadratures  rules  to  integrate  band-limited 
functions  is  roughly  7r/2  times  less  than  the  numbers  of  Gaussian  nodes;  the  numbers 
of  nodes  required  by  our  interpolation  formulae  in  order  to  interpolate  band-limited 
functions  is  roughly  7r/2  times  less  than  the  number  of  Chebychev  (or  Gaussian)  nodes. 
Again,  the  dependence  of  the  required  number  of  nodes  on  the  accuracy  requirements  is 
weak. 

4.  The  norm  of  the  differentiation  operator  based  on  our  nodes  is  of  the  order  as 
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compared  to  the  norm  of  the  spectral  dijfferentiation  operators  obtained  from  classical 
polynomial  expansions;  this  might  be  useful  in  the  design  of  spectral  (or  pseudospectral) 
techniques. 


9  Miscellaneous  Properties 

Prolate  spheroidal  wave  functions  possess  a  rich  set  of  properties,  vaguely  resembling  the 
properties  of  Bessel  functions.  This  section  establishes  some  of  those  properties.  Some 
of  the  identities  below  can  be  found  in  [20], [17], [5];  others  are  easily  derivable  from  the 
former. 


The  identity 

00 

^icxt  ^  ^  (239) 

j=o 

(see  Section  5)  has  a  number  of  consequences  which,  while  fairly  obvious,  seem  worth 
recording,  since  similar  properties  of  other  special  functions  have  often  been  found  useful. 
Differentiating  (139)  m  times  with  respect  to  x  and  n  times  with  respect  to  t  yields  the 
formula 

/  1  \  (m+n)  oo 

(-)  (140) 

J=0 


for  all  x,t  e  [-1,1]-  Multiplying  (139)  by  and  integrating  with  respect  to  t, 

converts  it  into 


sin{c  •  (x  —  u)) 
X  —  u 


r  °° 
^  i=0 


(141) 


Taking  the  squared  norm  of  (139),  and  integrating  with  respect  to  x  and  t,  yields  the 
formula 


i=o 

combining  this  with  (21)  yields 
i=o  ^ 


(142) 


(143) 
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Setting  x  =  t  =  l  converts  (139)  into 


00 


j=o 

(144) 

The  identity 

• 

Xjipj{x)  =  J  dt 

(145) 

(see  Section  2.5)  also  has  a  number  of  simple  but  potentially  useful  consequences.  Dif¬ 
ferentiating  it  k  times  with  respect  to  re,  we  get 

J  d.t. 

(146) 

We  next  consider  the  integral 

f{x)  =  fia,x)  =  /  ^  ^  _  ^^i(^)  dt. 

(147) 

Differentiating  (147)  with  respect  to  x,  we  have 

(148) 

Multiplying  (147)  by  ica,  and  subtracting  it  from  (148),  we  obtain 

■^f{a,x)  —  icaf{a,x)  =  ic  J  ipj{t)  dt 

(149) 

=  icXjipj{x). 

In  other  words,  /  satisfies  the  differential  equation 

f'{x)  —  icaf(x)  =  icXjipj{x). 

(150) 

The  standard  “variation  of  parameter”  calculation  provides  the  solution  to  (150): 

fix)  =  icXj  r  dt  +  /(O) 

J  0 

(151) 

Introducing  the  notation 

^  1  d 

ic  dx 

(152) 
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(i.e.  V  is  the  product  of  multiplication  by  l/ic  and  differentiation),  we  rewrite  (146)  as 

=  ^  (153) 

Aj  J-i 

for  an  arbitrary  polynomial  P  (with  real  or  complex  coefficients), 

P{V){^i)(x)  =  y  £  P(t)  (154) 

By  the  same  token,  the  function  4>  defined  by  the  formula 

satisfies  the  differential  equation 

P{V){<f>)ix)  =  XraXPmix).  (156) 

The  following  lemma  provides  a  recursion  connecting  the  values  of  the  k-th.  derivative 
of  the  function  'ipm  with  its  derivatives  of  orders  k  -1,  k  -2,  k  -Z,  k  -  4. 

Lemma  9.1  For  any  positive  real  c,  integer  m  >  0,  and  x  G  (— oo,  +oo), 


(1  -  x^)  -2{k^l)x  ^i^+^)(x) 

+  {Xm-k{k  +  l)-  ^  x^)  '4>^J^{x) 

-2(?kx  -c^  k{k-l)  ipl^~'^\x)  =  0  (157) 

for  all  k>  2.  Furthermore, 

(1  -  x^)  ^'^{x)  -  43:  V’m(^)  +  (Xm  -  2  -  x^)  xp'^{x) 

-  2  X 'ipmix)  =  0 .  (158) 

In  particular, 

-2ik  +  l)  V-r‘>(l)  +  (Xm  -*:(«:+  1)  -  <?)  i’S'll) 

-2^k i/-£-‘>(l)  -c^k(k-\)  =  0  (159) 

for  all  k  >2,  and 

-2xl;'^{l)  +  {Xm-c^)Ml)  =  0,  (160) 

-  4 1/;"  (1)  +  (x^  -  2  -  c^)  ^'^{1)  -  2  ^^(1)  =  0 .  (161) 
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Furthermore,  for  all  integer  m  >  Q  and  k  >2, 

-cH{k-l)i,i^-^\0)=0.  (162) 

For  all  odd  m, 

C(0)  +  (x™-2)V’;„(0)  =  0,  (163) 

and  for  all  even  m, 

^m(O)  +XmV^m(0)  =  0.  (164) 

Finally,  for  all  integer  m  >  0,  /c  >  0, 

li'™(l)#0,  (165) 

4“li(0)=0,  (166) 

./.g+''(0)  =  0.  (167) 

Proof.  All  of  the  identities  (157)  -  (164),  (166),  (167),  are  immediately  obtained  by 
repeated  differentiation  of  (24). 

In  order  to  prove  (165),  we  assume  that 

t^m(l)  =  0  (168) 

for  some  integer  m  >  0,  and  observe  that  the  combination  of  (168)  with  (159),  (160),  (161) 
implies  that 

V-S>(1)  =  0  (169) 

for  all  A:  =  0, 1, 2, ... .  Due  to  the  analyticity  of  iii  the  complex  plane,  this  would 
imply  that 

ipm{x)  =  0  (170) 


for  all  a:  G 

The  following  is  an  immediate  consequence  of  the  identity  (160)  of  Lemma  9.1. 
Corolleiry  9.2  For  all  integer  m,  n  >  0, 

V'm(l)  •  '*An(l)  -  ■0n(l)  *  V'm(l)  =  (Xn  "  Xm)  '  'l/’n(l)  '  V’m(l)  , 


□ 


(171) 


where  XmiXn  €  R  oJ'e  as  defined  in  Theorem  2.6. 


Theorem  3.1,  in  Section  3.1,  gives  formulae  for  the  entries  of  matrices  for  differen¬ 
tiation  of  prolate  series  and  for  multiplication  of  prolate  series  by  x.  Matrices  for  any 
combination  of  differentiation  and  of  multiplication  by  a  polynomial  can  obviously  be 
constructed  from  these  two  matrices;  for  instance,  calling  the  differentiation  matrix  D, 
and  the  multiplication- by-rr  matrix  X,  the  matrix  for  taking  the  second  derivative  of  a 
prolate  series,  then  multiplying  it  by  5  —  rc^,  is  equal  to  (5/  —  X^)D‘^. 

In  many  cases,  however,  there  are  simpler  formulae  for  the  entries  of  such  matrices, 
that  is,  for  inner  products  of  ‘ipj{x)  with  its  derivatives  and  with  polynomials.  The  follow¬ 
ing  theorems  establish  several  such  formulae,  as  well  as  a  few  formulae  for  inner  products 
which  do  not  involve  ipjix)  itself  but  only  its  derivatives.  We  start  with  Theorem  3.1, 
restated  here  for  consistency. 

Theorem  9.3  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm  =  n  (mod  2),  then 

j  ^  ■0n(^)  dx  =  j  ^  Xlpnix)  'Ipmix)  dx  =  0.  (172) 

If m^n  (mod 2),  then 

[  'tpnix)  Ipmix)  dx  =  ^m(l)^n(l),  (173) 

•/-I  +  K 

f  X-lpn{x)lpm{x)  dx  =  j-  ■^-^2  V>m(l)  V>n(l)  •  (174) 

J-i  -1-  XI 

Theorem  9.4  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod 2),  then 

J  ^xi;'^{x)'tlJm{x)  dx  =  0.  (175) 

If  m  =  n  (mod  2),  then 

[  X1p'^{x)'llJm{x)  dx=  (2V’m(l)^n(l)  “  <^mn)  •  (176) 

J —  1  Afn  I  An 

Proof.  Identity  (175)  is  obvious  since  the  functions  ipj  are  alternately  even  and  odd  (see 
Theorem  2.4).  In  order  to  prove  (176),  we  consider  the  integral 

J  ^xil}'^{x)'ipm{x)  dx 

=  Y  I  U  l  '^rn{x)  dx 
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=  ^  J  ^xrprnix)  (^j  dt^  dx 

^  ^  I-l  ^  (Xl  ^ 

J— 1 


In  other  words, 


fi  Am  r^ 

/  Xtp!^{x)i}rnix)dx=-^  X  ijnix)  dx  . 

On  the  other  hand,  integrating  the  left  side  of  (177)  by  parts,  we  obtain 
X  iIj'^{x) 'ipmix)  dx 

=  2'l/jm(l)llJn{l)  -  J_^  i'lpriix)  Ip'mix)  X  +  1pn{x)  1pm{x))  dx 
=  2ll;rn{l)‘tpn{l)  -  X  xl}rt{x)  1p'^(x)  dx  -  5mn  ■ 
Combining  (177)  and  (178),  we  have 

X'ip'^(x)'ijjn(x)  dx 
An 

=  2lprn{l)  ^n(l)  -  ^  X  ^m(^)  ^n(a^)  dx  -  6mn  , 

from  which  (176)  follows  directly. 


(177) 


Theorem  9.5  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod 2),  then 


J  ^  X^  1p'^{x)  Tpmix)  dx  =  0. 

If  m  =  n  (mod  2)  and  n,  then 

i-i 

J_  dx 


(178) 


Am  An 


An  “1“  An 


(V’n(l)^m(l)  -'i/’m(l)^n(l)) 

-'0n(l)  V’m(l) 


(179) 
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(180) 


Ati 


Am  An 
4A„ 


(Xn  Xm)  ^n(l)  ^m(l) 


l/»„(l)V'm(l), 


An  +  K 

where  Xm,Xn  ^  o,i"e  as  defined  in  Theorem  2.6. 


Proof.  Clearly  (178)  is  true,  since  the  functions  ipj  are  alternately  even  and  odd.  In 
order  to  prove  (179)  and  (180),  supposing  that  m  =  n  (mod  2)  and  m  ^  n,  we  consider 
the  integral 

rl 

'Ipmix)  dx 

=  Y  '  iI-1 

=  -Y  '>Pm(x)  X^  •  dx 

=  -Y  ^m(a:)rr^e*‘=^*dx^  ?/>n(t)t^  dt 

=  ^  [  t^V'n(t)  C(<)  dt, 

Aji  J-1 

which  is  summarized  as 

'Ip'^ix)  Iprnix)  dx  =  ^  X^  1p'^{x)  'Ipnix)  dx  .  (I8l) 

On  the  other  hand,  integrating  the  left  side  of  (181)  by  parts,  we  have 
X^  1p'^{x)  Iprnix)  dx 

=  2V'n(l)V'm(l)  i’nix)  (i^'^ix)  X^  +  2  X 'lpm{x))  dx 

=  2  V'n(l)  V'm(l)  “  2  V’n(^)  X  dx 

-  dx.  (182) 

Due  to  Theorem  9.4  and  the  fact  that  m^n,we  immediately  rewrite  (182)  as 

rl 

1^  X^  1p'^{x)  1pm{x)  dx 

2A 

=  2^;(l)^m(l)-  .  '  2^n(l)V^m(l) 

'^m  ■> 

-  dx,  (183) 
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which  we  rewrite  as 

dx 

=  2l/>;,(l)V>m(l)  -  T^-^r  V>n(l)^m(l) 

-Ajtx  “T  Aji 

-  j_^x^il}l{x)il)m{x)  dx.  (184) 

Swapping  m  with  n,  we  convert  (184)  into 

'^nix)  dx 

4  A 

=  2C(l)V^n(l)  -  T  V>n(l)^m(l) 

I 

-  j  ^  ip'^(x)  ipnix)  dx.  (185) 

Combining  (184)  and  (185),  we  obtain 

[  x'^  ^"(x)  ^rnix)  dx-2  Xp!^{l)  V’m(l)  + 

>/  — 1  Ajji  “T  Aji 

=  [  z^V»"(a:)V'n(2;)  c?a:-2^(„(l)^„(l)+ — ^^^^„(l)^„j(l),  (186) 

•'—1  A,tj  +  A,j 

which  is  obviously  equivalent  to 

x"^  ipn{x) 'tpmix)  dx 

=  V'n(2:)  da;  +  2  (^’(,(1)  l/’m(l)  -  V’m(l)  V’n(l)) 

+ 4  y"  , 'r  • 

Aji  ~r  Atxi 

Finally,  combining  (181)  with  (187),  we  have 

^  f  (x)^„(x)  dx 

An  J—l 

=  j_^  i’nix)  dx  +  2  (^;(1)  i)ra{l)  “  l/'m(l)  V’n(l)) 

+  4  ^  ^  V’n(l)  ^m(l)  , 

An  “T  Ajn 


(187) 


which  is  easily  rewritten  as 


=  2«(1)^™(1)-V’;.(l)^n(l)) 


A. 


+  4 


Xfi  At] 

An  +  A;] 


or 


2  An 


Am  An 
4  A„ 


An  +  At 


«(l)V’m(l)-C(l)^7.(l)) 

V'n(l)^m(l)- 


(188) 

□ 


We  finally  rewrite  (188)  as  (180)  using  Corollary  9.2. 

The  following  theorem  is  an  immediate  consequence  of  combining  the  preceding  theorem 
with  equation  (184)  from  its  proof. 

Theorem  9.6  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod 2),  then 


(te  =  0. 


(189) 


•  Ifm  =  n  (mod  2)  and  m^n, 

O  \ 

=  2^;(1)^„(1)  +  3— (V.;»(1)V’„(1)-C(1)W1))  (190) 

^  Am  An 

O  \ 

=  2*(l)V’„(l)  +  y-f^(*(l)V>„(l)-C(l)V'«(l))  (191) 

=  ^„(i)^„(i)  I  (192) 

Theorem  9.7  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
-  non-negative.  Ifm^n  (mod 2),  then 

Li  dx  =  j  ^X^  Xpnix)  ^m{x)  dx  =  0  (193) 
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Ifm  =  n  (mod  2)  and  mj^n,  then 

=  3;r^«(i)'/’».(i)-'A„(i)!^;(i))  (194) 

=  (Xn  -  Xm)  V’m(l)  <^.(1),  (195) 

j  ^X^1pri{x)'ipm{x)  dx 

=  (186) 

=  ~4  \r~\2  ~ 

-  K 

where  XmjXn  €  o,re  as  defined  in  Theorem  2.6. 

Proof.  Identity  (193)  is  obvious,  since  the  functions  xpj  are  alternately  even  and  odd. 
In  order  to  prove  (194)-(197),  we  start  with  the  expression 

Xnipnix)  =  -c^  J  dt.  (198) 


Taking  the  inner  product  of  (198)  with  ipmix),  we  have 
An  J  ^  xp'^ix)  Xpm{x)  dX 

=  J  ^  (^J  ^  dtj  ipm{x)  dx 

=  J  ^  (^J  ^  1pm{x)e^‘^''^  dx'^  dt 

fl 

=  -C^Xm  J  ^  dt, 

which  we  summarize  as 

x^  i)n{x)  ■0r7i(a:)  dx  =  ~  <(x)  'ipmix)  dx .  (199) 

Swapping  n,  m,  we  rewrite  (199)  in  the  form  of 

j  ^X^'llJn{x)'ll^m{x)dx 

=  “4^  /  ‘^mix)  tpn{x)  dx .  (200) 

An  J-1 
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(201) 


Combining  (199)  and  (200),  we  get 

^m(^)  i’nix)  dx  . 

On  the  other  hand,  integrating  the  left  side  of  (201)  by  parts,  we  have 
j_^^n{x) 'tpmioc)  dx 

=  i’M  ^m(l)  -  ^^,(-1)  1pm{-l)  -  'tpnix)  C(^)  dx 

=  2  t/);(l)  i^rnil)  -  (^n(l)  C(l)  “  V’n(-l)  ^m(-l)) 

+  'ipnix)  ip'^{x)  dx.  (202) 

We  rewrite  (202)  in  the  form  of 
J_^'>Pnix)  iptnix)  dx 

=  2  (^’^(1)  V>m(l)  -  '4>n{l)  C(l))  +  J_^  M^)  C(^)  dx  . 

We  combine  (201)  and  (203)  and  get 

- 1)  /j  w  cw  * 

=  ■  (203) 

Since  m  ^  n,  we  easily  rewrite  (203)  as  (194).  We  obtain  expression  (196)  by  combin¬ 
ing  (200)  and  (194).  The  identities  (195),  (197)  follow  from  (194),  (196)  immediately 
due  to  Corollary  9.2.  □ 


Theorem  9.8  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
%  non-negative.  Let 


ry 

^n(y)  =  /  i’nix)  dx. 

Jo 

(204) 

If  n  is  odd  and  m  is  even,  then 

1 

• 

J  ^  -  ■ijjnit)  1pm(t)  dt 

(205) 

- 

(206) 

(207) 
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(208) 


Ifm  —  n  (mod  2),  then 

j'tpn{t)rpm(t)  dt  =  0. 


Proof.  We  start  with  the  identity 

Ktpnix)  =  j  dt. 

Integrating  (209)  with  respect  to  x,  we  have 

ry 

An  /  ^n{x)  dx 

Jo 

(209) 

=  dtj  dx 

(210) 

=  J  'ipni't)  j  dxdt 

(211) 

which  we  summarize  as 

(212) 

An  ^n{y)  =  -  /'  7  Mt)  dt-^  f'  1  ^n(t)  dt . 

IC J-l  t  tc J-1  t 

Taking  the  inner  product  of  (213)  and  ipmiy),  we  obtain 

(213) 

K  J ^'^n{y)i^m{y)  dy 
~  Tcl  1  '{Ill 

(214) 

=  h  /-1 1  ■  (/-i  ‘'®')  * 

1  ri  1  /-i 

-—  /  7V’n(i)  c?t-  /  V’m(2/)  dy 

%  C  J  — 1£  »/ — 1 

(215) 

=  7V'n(i)^m(t)  dt 

^C  J-l  t 

• 

1  /•!  1  /•! 

-—  /  7^n(^)  dt  -  /  ^„j(y)  dy, 

2CJ-1  t  J-l 

(216) 
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which  we  summarize  as 


=  /  ^n{t)^lJm{t)  dt 

•'“1 

+  ^  f  ];^n{t)dt-f  'ipmiy)  dy  (217) 

*'“’1  ^  */“l 

Exchanging  m  with  n,  we  convert  (217)  into 

j  'tpn{t)  dt 

=  ic^  [  ^rn{t)'4’n{t)  dt 
An  J 

+  YJ_^-i}mit)dt‘j^^tpniy)dy,  (218) 

and  combining  (217),  (218),  we  get 

^ic  f  '^n{t)  V’m(i)  dt-^ic  [  -^rnit)  i>n{t)  dt 

Am  •'—1  An  •'—1 

=  -p  /  7^m(^)  dt-  f  Xpn{t)  dt 
An  J-l  t  J-1 

f  ^i^n{t)dt-f  'lprnit)dt.  (219) 

Ayyj  «/— 1  T  J ~“X 

Suppose  that  m  is  even  and  n  is  odd;  then  the  first  product  in  the  right  hand  side  of  (219) 
is  zero,  so 

^ic  [  ^n(i)  ■tpmit)  dt-^ic  f  i>n{t)  dt 

Am  -/-I  A„  J-l 

=  /  ^lpnit)dt-f  ‘lj;m{t)dt,  (220) 

*'“1  t  J '-1 

which  is  equivalent  to 

j^^'i!n{t)'lpm{t)  dt 

=  %■  /  '^m{t)lpn{t)  dt 

A„  -/-I 

1  1  /•!  1  /•! 

-  —  —  /  -  1pn{t)  dt  ■  /  Iprnit)  dt , 

A72  %  C  J “-X  Tf  «/  — 1 


(221) 


or 


=  dt 

+  7 '>Pnit)  dt  •  f  1prn{t)dt.  (222) 

^  C  •/ —1  X  v  “1 

On  the  other  hand,  integrating  the  left  side  of  (222)  by  parts,  we  obtain 
J^^'^rnit)'lpnit)  dt 

=  ^„(1)  ^r^(l)  -  ^„(-l)  ^n(i)  dt  .  (223) 

Since  the  product  ^m(^)  ^71(2;)  is  an  odd  function  when  m^n  (mod  2),  we  rewrite  (223) 
as 

f  ^m(t)  V’n(i)  dt 

=  2^  n{l)^  mil)- f_^'^n{t)Mt)dt.  (224) 

The  combination  of  (222)  and  (224)  implies  that 

f  ^nit)  'Ipmit)  dt+^  [  ^nit)  'tpm{t)  dt 

=  2  ^nil)  ’®^m(l)  f  7  V’n(t)  dt-f  ^rn{t)  dt  ,  (225) 

^rn  XCJ—1  t  J—1 

or 

\2  I  \2  -1 

<lfnit)^Pmit)dt 

•'-1 

=  2^^„(l)^f^(l)  7^n(i)  dt-  f  iprnit)  dt ,  (226) 

icJ-i  t  J-l 

which  is  equivalent  to 

J_^^nit)'(pm{t)  dt 
0  \2 

" a!^  h  L  \  W  *  •  /-. 
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M 


Finally,  combining  (217)  and  (227),  we  have 
ri  1 


2  A,7i  A„ 


=  ic 


A^  +  A^ 
Am 

■A2+A2 


^n(l)^m(l) 


^  tpnit)  dt  ■  ^m(t)  dt . 


(228) 


Equation  (208)  is  easily  proven  since  the  product  |  '^n{x)  is  an  odd  function  when¬ 
ever  m  =  n  (mod  2).  □ 

The  above  theorems  do  not  use  much  of  the  detailed  structure  of  the  integral  operators 
of  which  the  functions  {'ijjj}  are  eigenfunctions.  Thus  many  of  them  generalize  easily  to 
the  case  of  an  operator  L  :  L'^[0, 1]  -4  ^^[0, 1]  defined  via  the  formula 

L{v){x)  =  f  K{xt)-ip(t)  dt,  (229) 

Jo 

for  some  function  K  :  [0, 1]  ->  C;  the  following  theorem  is  an  example  of  this. 

Theorem  9.9  Let  Ai.Aj  be  two  eigenvalues  of  the  operator  L  defined  by  (229),  that  is, 


Then 


[  K {xt) ‘tl)i{t)  dt  =  AiV'i(a;), 
Jo 

[  K {xt) -11)2(1)  dt  =  X2'tp2ix). 
Jo 

Aj  xi;[{x)-ip2{x)  dx 


(230) 

(231) 

(232) 


[  X -if'^ix) 'tpi{x)  dx 
Jo 

provided  that  neither  Ai  nor  the  denominator  of  the  right  hand  side  of  (232)  is  zero. 
Proof.  Differentiating  (230),  (231)  with  respect  to  x,  we  get 

[  tK'{xt)'il)i{t)  dt  —  Xiif[{x),  (233) 

Jo 

[  tK'{xt)'i/)2{t)  dt  =  X2'4>2{x).  (234) 

Jo 

Multiplying  (233)  by  X'if2{x),  we  have 

XiX'il)[{x)'if2{x)  =  x'il)2{x)  [  tK'{xt)ip\{t)  dt.  (235) 

Jo 
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Integrating  on  the  interval  [0, 1],  we  obtain 


[  x'ip[{x)'ilj2{x) 
Jo 

dx  = 

J  X'ip2ix)  J  tK'{xt)ipi{t)  dtdx 

(236) 

= 

J  i'fpiit)  J  X  K' {xt)  -02  (a:)  dx  dt. 

(237) 

Renaming  the  variables  of  integration  Or 

i  the  right  hand  side  from  xtot  and  vice 

versa. 

we  get 

-^1 

[  XJp[{x)'ip2ix) 

fO 

dx  = 

J  x'(pi{x)  j  t K' (xt) 'ip2it)  dtdx. 

(238) 

Substituting  (234)  into  (238),  we 

obtain 

[  Xil;[{x)^2ix) 
lo 

dx  = 

X2  [  x'ipi{x)7p2{x)  dx, 

Jo 

(239) 

from  which  (232)  follows  immediately,  as  does  its  caveat.  □ 


The  following  theorem  establishes  the  relation  between  the  norm  of  each  function  ■ipj 
on  [—1, 1]  (which  in  this  paper  is  taken  to  be  one),  and  its  norm  on  (— oo,  oo). 

Theorem  9.10  Suppose  that  c  is  real  and  positive,  and  that  the  integer  n  is  non¬ 
negative.  Then 


/OO  1 

dx=  — 

-oo  Pn 


(240) 


where  is  given  by  (21). 

Proof. 

roo 


/oo 

-00 


dx  = 


J —oo  \7rfJ,n  1—i 

=  -[\n{t)-(-r 

P'n  \  TT  J- 

= 

Un 


X  —  t 

1  f°°  sin(c  •  (x  —  t)) 

X  —  t 


dtj  dx 

ipn{^)  da:']  dt 


P'n 

A^n 


□ 


The  following  theorem  extends  Theorem  (9.10)  to  any  band-limited  function  with 
band  limit  c. 
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Theorem  9.11  Suppose  that  c  is  real  and  positive,  that  the  integer  n  is  non-negative, 
and  that  f  :  R <D  is  a  band-limited  function  with  band  limit  c.  Then 


Proof. 


/oo  1 

ifnix)  f{x)  dx  =  —  /  f{x)  dx. 

-OO  {Jin  •'-1 


/oo 

Ipnix)  f{x)  dx 

-oo 

=  [  (—  [  -  fix)  dx 

J-oo  y-irpn  J-l  X  —  t  J 

=  —  /  i’nit)  •  (-  f  (if 

fJiji  */  — 1  ^TT  J—oo 

=  —  f  i^ni't)fit)  dt. 

Un 


X  —  t 


(241) 


□ 


Theorem  9.12  Suppose  that  c  is  real  and  positive,  and  that  the  integer  n  is  non¬ 
negative.  Then 


(OO  I  —i’mix),  if  -  I  <  X  <  1, 

f  =  J  ^ 

J  — OO 

0,  if  X  >  1  or  X  <  —1. 


(242) 


Proof.  Since  xfrn  is  an  eigenfunction  of  the  operator  Qc  defined  in  (19),  and  fj,m  is  the 
corresponding  eigenvalue, 


1 

/^m^m(^)  “  / 

TT  J-1 


1  /•!  sm(c  •  (x  —  u)) 


x  —  u 


du. 


(243) 


Thus 


/OO 

e'‘^‘Mt)dt 

-OO 

1  gicrt  (^1  sin(c-(x-u)) 
pm  J —'X 


X  —  U 


J-  ['  ( i  r  e-  *)  dn 

Pm  J-i  V’’’  X  —  U  J 


dt 

(244) 

dw 

(245) 
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Since  the  innermost  integral  is  the  orthogonal  projection  operator  onto  the  space  of 
functions  of  band  limit  c  on  (— oo,cx)),  applied  to  the  function  it  follows  that; 


/OO 

-OO 


1 

=  -  /  V’m(M) 

Um  •'“1 


if  -  1  <  X  <  1, 

0,  if  a;  >  1  or  a:  <  —1 


1  /•! 


J  i’mi'U')  du,  if  -  1  <  a:  <  1, 


(246) 


(247) 


if  a:  >  1  or  x  <  — 1, 


from  which  (242)  follows  immediately. 


The  following  five  theorems  establish  formulae  for  the  derivatives  of  prolate  functions 
and  their  associated  eigenvalues  with  respect  to  c. 


Theorem  9.13  For  all  positive  real  c  and  non-negative  integer  m, 

2t^^(l)  -1 

dc  2c  ' 

Proof.  We  start  with 


AmV'm(a:)  =  'Iprnit)  dt . 


(248) 


(249) 


Differentiating  (249)  with  respect  to  c,  we  obtain 
dXm  ,  r  ^  y  dxpm{x) 


=  J  dt  +  j  e*® 


dlpm{t) 


(250) 


Multiplying  by  on  both  sides  of  (250),  and  integrating  on  the  interval  [—1, 1],  we 

get 


=  J  y  dt  dx 


(251) 


j  ^  ipm{x)  J  ^  dt  dx . 
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which  we  rewrite  as 


d^m  ,  . 


L 


dc 


'lpm{x)  dx 


=  J  itlpm{t)  j  X -Iprai^)  dx  dt 

d'lpmjt)  ^icxt . 

J-1  dc  J-i 


-ipnix)  dx  dt 


—  I  ^  ^  '0m  (^)  .  ^  dt 

J-1  ic  dt 


+  A, 

which  we  summarize  as 
dXm 


I 


dc 


Ipriiit)  dt, 


dc 


Am  j  dlpfjiit) 


On  the  other  hand,  integrating  the  right-hand  side  of  (254)  by  parts,  we  have 


=  0m(l)  +  0m(-l)-l-/j0m(i)i 


dt 


d'lpmjt) 

dt 


dt . 


which  we  rewrite  as 


^  dt  =  V'J.(l)  -  i 


Finally,  substituting  (256)  into  (254),  we  get 


dXrrt 

dc 


A. 


20^(1) -1 

2c 


(252) 


(253) 


(254) 


(255) 


(256) 

(257) 

□ 


Theorem  9.14  For  any  positive  real  c  and  non-negative  integer  m, 
dpm  2  i2 


(258) 
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Proof.  We  start  with  the  identity 


-2ct  , 

A^m  —  ^rfL  '^m- 
TT 


DiflFerentiating  (259)  with  respect  to  c,  we  get 


9fi 


2c  f  ^  dXr 


m  _  I  —  '-"'m  ,  \ 

““  I  r\  I 


dc 


dc  TT 

Substituting  Lemma  9.13  into  (260),  we  get 
dll. 


"b  Xm  Xr 


TT 


dc 


?£.2a  a  ^JMz1  +  1x  X 

TT  2C  TT 


^  2t/;2(l)-l  1 

—  ^  llm - 7: - 1 - llm 

2c  C 

2  /2  /.X  1  1 

—  IlmWmK^j  jlm  llm 

C  C  C 

2 

=  -A^m^^(l). 

O 


(259) 


(260) 


(261) 


(262) 

□ 


The  following  theorem  immediately  follows  from  Theorems  9.13  and  9.14. 

Theorem  9.15  For  all  positive  real  c  and  non-negative  integer  m,  n, 

=  ■  (264) 

\  l^n  /  l^n  ^ 

Theorem  9.16  Suppose  that  c  is  real  and  positive,  and  the  integers  m,n  are  non¬ 
negative.  If  m  ^n,  then 

^^•0m(i)^(i)  dt  = -^^^^^3^V’m(l)V'n(l)-  (265) 

If  m  =  n,  then 

l\m{t)^{t)  dt  =  0.  (266) 
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Proof.  Since  the  norm  of  i()n  on  [—1, 1]  remains  constant  as  c  varies,  must  be  orthog¬ 
onal  on  [—1, 1]  to  its  own  derivative  with  respect  to  c,  which  immediately  yields  (266). 
To  establish  (265),  we  start  with  the  identity 

K‘fpn{x)=[  ipriit)  dt .  (267) 

DiflFerentiating  (267)  with  respect  to  c,  we  get 


+  -^n 


d-lpn 

dc 


=  X  t  Mt)  +  ■  (268) 

Multiplying  both  sides  of  (268)  by  iptnix)  and  integrating  with  respect  to  x,  we  have 
An  'Ipmix)  dx 

=  ^  y  ^  X  i)^{x)  dx-\-Xrn  j  ^  V’m(i)  “ dt ,  (269) 

which,  using  (176),  we  rewrite  as 


(A„  -  Xm)  dt 


Xfi  A^ 


(2^„j(l)V’n(l)  ^mn) 


C  An;  +  An 

Assuming  that  m  ^  n,  and  dividing  by  An  -  Am,  we  then  get  (265). 


(270) 

□ 


Theorem  9.17  Suppose  that  c  is  real  and  positive,  and  the  integer  m  is  non-negative. 
Then 

^  =  2c_y^x2^^(x).  (271) 

Proof.  Due  to  Theorem  2.6, 

(1  -  X^)V’"  (x)  -  2x^(„(x)  +  iXm  -  (?X^)  1prn{x)  =  0.  (272) 

Making  the  infinitesimal  changes  c  =  c  h,  Xm  =  Xm  +  s,  and  'tprn{x)  =  '0m(x)  -\-  5(x), 
this  becomes 

(1  -  x^)  •  (V>"  (x)  -h  5"(x))  -  2x  ■  {ip'mix)  H-  5'(x)) 

+  (Xm  +  £  -  (c  +  h)  V)  ■  {'ipmix)  +  <J(x))  =  0.  (273) 
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Expanding  each  term,  discarding  infinitesimals  of  the  second  order  or  greater  (that  is, 
products  of  two  or  more  of  the  quantities  7i,  £,  and  6{x)),  and  subtracting  (272),  we  get 


(1  —  x^)  S"(x)  —  2xS'(x)  4-  (Xm  -  +  (£■  ~  2chx^)'^rn{x)  =  0.  (274) 

Let  the  self-adjoint  differential  operator  L  be  defined  by  the  formula 

L{f){x)  =  {l-x^)f"{x)-2xf{x)-\-{xm-c^x^)f{x).  (275) 

Then,  multiplying  (274)  by  'il;m(x)/h  and  integrating  on  [-1, 1],  we  get 

Li  Li 

Now  I  =  In  addition,  since  L  is  self-adjoint, 

/ 1  dx  =  HMix)  dx.  (277) 

But  due  to  (272),  L{xl'm){x)  =  0  for  all  x  €  [—1,1],  so  the  integral  (277)  is  zero; 
Thus  (276)  becomes 

^  =  2c/_‘ x\-i(x).  (278) 

□ 


10  Generalizations  and  Conclusions 

In  this  paper,  we  design  quadrature  rules  for  band-limited  functions,  based  on  the  prop¬ 
erties  of  Prolate  Spheroidal  Wave  Functions  (PSWFs),  and  the  connections  of  the  latter 
with  certain  fundamental  integral  operators  (see  (17),  (19)  in  Section  2.5).  The  quadra¬ 
tures  are  a  surprisingly  close  analogue  for  band-limited  functions  of  Gaussian  quadratures 
for  polynomials,  in  that  they  have  positive  weights,  are  optimal  in  the  appropriately  de¬ 
fined  sense,  and  their  nodes,  when  used  for  approximation  (as  opposed  to  integration), 
result  in  extremely  efficient  interpolation  formulae.  Thus,  Sections  5-7  of  this  paper  can 
be  viewed  as  reproducing  for  band-limited  functions  much  of  the  standard  polynomial- 
based  approximation  theory  (for  which  see,  for  example,  [24]).  Generally,  there  is  a 
striking  analogy  between  the  band-limited  functions  and  polynomials. 

Obviously,  there  are  certain  differences  between  the  resulting  apparatus  and  the  stan¬ 
dard  numerical  analysis.  To  start  with,  where  the  classical  techniques  are  optimal  for 
polynomials,  the  approach  of  this  paper  is  optimal  for  band-limited  functions;  whenever 
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the  functions  to  be  dealt  with  are  naturally  represented  by  trigonometric  expansions  on 
finite  intervals,  our  quadrature  and  interpolation  formulae  tend  to  be  more  efficient  than 
those  based  on  the  polynomials.  When  the  functions  to  be  dealt  with  are  naturally  rep¬ 
resented  by  polynomials,  the  classical  approach  is  more  efficient;  however,  many  physical 
phenomena  involve  bend-limited  functions,  and  very  few  involve  polynomials. 

Qualitatively,  the  quadrature  (and  interpolation)  nodes  obtained  in  this  paper  behave 
like  a  compromise  between  the  Gaussian  nodes  and  the  equispaced  ones:  near  the  middle 
of  the  interval,  they  are  very  nearly  equispaced,  and  near  the  ends,  they  concentrate 
somewhat,  but  much  less  than  the  Gaussian  (or  Chebychev)  nodes  do.  For  large  c,  the 
distance  between  nodes  near  the  ends  of  the  interval  is  of  the  order  with  the  total 
number  of  nodes  close  to  -.  In  contrast,  the  distance  between  the  Gaussian  nodes  near 
the  ends  of  the  interval  is  of  the  order  with  n  the  total  number  of  nodes.  A  closely 
related  phenomenon  is  the  reduced  norm  of  the  differentiation  operator  based  on  the 
prolate  expansions:  for  an  n— point  differentiation  formula,  the  norm  is  of  the  order 
as  opposed  to  for  polynomial-based  spectral  differentiation.  Thus,  PSWFs  are  likely 
to  be  a  better  tool  for  the  design  of  spectral  and  pseudo-spectral  techniques  than  the 
orthogonal  polynomials  and  related  functions. 

Much  of  the  analytical  apparatus  we  use  was  developed  more  than  30  years  ago 
(see  [20]-[21],  [17],  [18]);  the  fundamental  importance  of  these  results  in  certain  areas  of 
electrical  engineering  and  physics  has  also  been  understood  for  a  long  time.  However, 
there  appears  to  have  been  no  prior  attempt  made  to  view  band-limited  functions  as  a 
source  of  numerical  algorithms.  Generally,  there  is  a  fairly  limited  amount  of  information 
in  the  literature  about  the  PSWFs,  especially  when  compared  to  the  wealth  of  facts  on 
many  other  special  functions.  Section  9  of  this  paper  is  an  attempt  to  remedy  this 
situation  to  a  small  degree. 

The  apparatus  built  in  this  paper  is  a  strictly  one-dimensional  one.  Obviously,  one 
can  construct  discretizations  of  rectangles,  cubes,  etc.  by  using  direct  products  of  one¬ 
dimensional  grids;  the  resulting  numerical  algorithms  are  satisfactory  but  not  optimal. 
Furthermore,  representation  of  band-limited  functions  on  regions  in  higher  dimensions 
is  of  both  theoretical  and  engineering  interest.  Obvious  applications  include  seismic 
data  collection  and  processing,  antenna  theory,  NMR  imaging,  and  many  others.  When 
the  region  of  interest  is  a  sphere,  most  of  the  necessary  analytical  apparatus  can  be 
found  in  [21].  At  the  present  time,  we  have  constructed  and  implemented  somewhat 
rudimentary  versions  of  the  relevant  numerical  algorithms;  we  are  conducting  numerical 
experiments  with  these,  and  will  report  the  results  at  a  later  date.  A  much  more  difficult 
set  of  questions  is  presented  by  the  structure  of  band-limited  functions  on  more  general 
regions. 
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Table  1:  Quadrature  performance  for  varying  band  limits,  for  e  =  10 


-7 


c 

n 

Maximum  Errors 

iVpol 

Roots 

Refined 

10.0 

9 

0.96E-05 

0.51E-07 

13 

20.0 

13 

0.17E-04 

0.94E-07 

19 

30.0 

17 

0.12E-04 

0.50E-07 

25 

40.0 

20 

0.70E-05 

0.30E-06 

31 

50.0 

24 

0.35E-05 

0.83E-07 

37 

60.0 

27 

0.25E-04 

0.27E-06 

43 

70.0 

31 

O.llE-04 

0.66E-07 

48 

80.0 

34 

0.48E-05 

0.17E-06 

54 

90.0 

38 

0.21E-05 

0.40E-07 

.  59 

100.0 

41 

0.12E-04 

0.91E-07 

65 

200.0 

74 

0.24E-05 

0.86E-07 

118 

300.0 

106 

0.32E-05 

0.21E-06 

171 

400.0 

139 

0.52E-05 

0.62E-07 

223 

500.0 

171 

0.56E-05 

0.88E-07 

275 

600.0 

203 

0.58E-05 

O.llE-06 

326 

700.0 

235 

0.57E-05 

0.12E-06 

377 

800.0 

267 

0.55E-05 

0.13E-06 

428 

900.0 

299 

0.53E-05 

0.14E-06 

479 

1000.0 

331 

0.50E-05 

0.14E-06 

530 

1200.0 

395 

0.44E-05 

0.13E-06 

632 

1400.0 

459 

0.38E-05 

O.llE-06 

734 

1600.0 

523 

0.31E-05 

0.97E-07 

835 

1800.0 

587 

0.28E-05 

0.80E-07 

937 

2000.0 

651 

0.23E-05 

0.64E-07 

1038 

2400.0 

778 

0.29E-05 

0.15E-06 

1240 

2800.0 

906 

0.19E-05 

0.84E-07 

1442 

4000.0 

1288 

0.37E-05 

0.17E-06 

2047 

54 


Table  2:  Quadrature  performance  for  varying  precisions,  for  c  =  50 


e 

Maximum  Errors 

Npoi 

Roots 

Refined 

O.lOE-01 

19 

0.45E-01 

O.lOE-01 

30 

O.lOE-02 

20 

0.70E-02 

0.13E-02 

32 

O.lOE-03 

21 

0.91E-03 

0.14E-03 

33 

O.lOE-04 

22 

0.82E-04 

0.13E-04 

34 

O.lOE-05 

23 

0.54E-04 

O.llE-05 

36 

O.lOE-06 

24 

0.35E-05 

0.83E-07 

37 

O.lOE-07 

25 

0.33E-05 

0.57E-08 

38 

O.lOE-08 

26 

0.18E-06 

0.36E-09 

39 

O.lOE-09 

26 

0.18E-06 

0.36E-09 

40 

O.lOE-10 

27 

0.17E-06 

0.21E-10 

42 

O.lOE-11 

28 

0.79E-08 

O.llE-11 

43 

O.lOE-12 

29 

0.78E-08 

0.56E-13 

45 

O.lOE-13 

30 

0.31E-09 

0.27E-14 

55 

55 


Table  3:  Interpolation  performance  for  varying  band  limits,  for  £  =  10  ^ 


c 

n 

Maximum  Errors 

A^poi 

Roots 

Refined 

Cheb. 

Leg. 

5.0 

13 

0.12E-06 

0.12E-06 

17 

17 

10.0 

18 

0.12E-06 

0.13E-06 

24 

25 

15.0 

22 

0.24E-06 

0.25E-06 

31 

32 

20.0 

26 

0.26E-06 

0.28E-06 

37 

39 

25.0 

30 

0.22E-06 

0.23E-06 

43 

45 

30.0 

33 

0.67E-06 

0.73E-06 

49 

51 

35.0 

37 

0.42E-06 

0.46E-06 

55 

57 

40.0 

41 

0.25E-06 

0.27E-06 

61 

63 

45.0 

44 

0.54E-06 

0.60E-06 

67 

69 

50.0 

48 

0.29E-06 

0.33E-06 

73 

75 

100.0 

82 

0.39E-06 

0.46E-06 

128 

131 

150.0 

115 

0.52E-06 

0.64E-06 

182 

186 

200.0 

147 

0.12E-05 

0.15E-05 

235 

239 

250.0 

180 

0.83E-06 

O.llE-05 

287 

292 

300.0 

212 

0.13E-05 

0.17E-05 

340 

345 

350.0 

245 

0.75E-06 

O.lOE-05 

392 

398 

400.0 

277 

O.lOE-05 

0.14E-05 

443 

450 

450.0 

309 

0.13E-05 

0.18E-05 

495 

502 

500.0 

341 

0.16E-05 

0.22E-05 

547 

554 

1000.0 

662 

0.16E-05 

0.24E-05 

1058 

1068 

1500.0 

982 

0.15E-05 

0.25E-05 

1566 

1578 

2000.0 

1301 

0.20E-05 

0.35E-05 

2072 

2086 

56 
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Table  4;  Interpolation  performance  for  varying  precisions,  for  c  =  25 


E 

n 

Maximum  Errors 

■Npoi 

Roots 

Refined 

Cheb. 

Leg. 

O.lOE-01 

21 

0.38E-01 

0.43E-01 

31 

34 

O.lOE-02 

23 

0.37E-02 

0.41E-02 

34 

36 

O.lOE-03 

25 

0.29E-03 

0.31E-03 

37 

39 

O.lOE-04 

26 

0.74E-04 

0.81E-04 

39 

41 

O.lOE-05 

28 

0.44E-05 

0.47E-05 

41 

43 

O.lOE-06 

30 

0.22E-06 

0.23E-06 

43 

45 

O.lOE-07 

31 

0.46E-07 

0.49E-07 

45 

47 

O.lOE-08 

32 

0.95E-08 

O.lOE-07 

47 

49 

O.lOE-09 

34 

0.36E-09 

0.38E-09 

49 

51 

O.lOE-10 

35 

0.67E-10 

0.70E-10 

51 

52 

O.lOE-11 

37 

0.21E-11 

0.22E-11 

53 

54 

O.lOE-12 

38 

0.36E-12 

0.37E-12 

54 

56 

O.lOE-13 

39 

0.59E-13 

0.63E-13 

98 

61 
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Table  5:  Quadrature  nodes  for  band-limited  functions,  with  c  =  50  and  e  —  10~^ 

This  table  contains  only  half  of  the  nodes  and  weights,  in  particular  those  for  which  the 
node  is  less  than  or  equal  to  zero;  reflecting  these  nodes  around  zero  yields  the  remaining 
nodes,  the  weight  for  the  node  at  —x  being  the  same  as  the  weight  for  the  node  at  x. 


Node 

-  .9904522459960804E-1-00 
-.9525601106643832E-1-00 
-.8927960861459153E-h00 
-.8186117530609125E+00 

-  .7350624131965875E-hOO 

-  .6452878027260844E-h00 

-  .5512554698695428E-1-00 
-.4542505281525226E-I-00 
-.3551568458127944E-h00 
-.2546173463813596E-h00 
-.1531287781860989E-t-00 
-.5110121484050418E-01 


Weight 

0.2413064234922188E-01 

0.5024347217095568E-01 

0.6801787677830858E-01 

0.7952155999100788E-01 

0.8706680708376023E-01 

0.9216240765763570E-01 

0.9569254015486106E-01 

0.9817257766311556E-01 

0.9990914516102242E-01 

0.1010880172648715E+00 

0.1018214308931439E+00 

0.1021735189986602E-f-00 
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Figure  2:  Maximum  norms  of  derivatives  of  prolate  spheroidal  wave  functions  for  c  =  200, 
and  of  normalized  Legendre  polynomials 


Table  6:  Quadrature  nodes  for  band-limited  functions,  with  c  =  150  and  e  =  10“^^ 

This  table  contains  only  half  of  the  nodes  and  weights,  in  particular  those  for  which  the 
node  is  less  than  or  equal  to  zero;  reflecting  these  nodes  around  zero  yields  the  remaining 
nodes,  the  weight  for  the  node  at  —x  being  the  same  as  the  weight  for  the  node  at  x. 


Node 


-.9982883010959975E-h00 

-.9911354691596528E-1-00 

-.9788315280982487E-I-00 

-.9621348937901911E-hOO 

-  .9418386698454396E+00 
-.9186509576802944E-1-00 

-  .8931541850293142E-1-00 
-.8658083894041821E-I-00 

-  .8369709588254746E-F00 
-.8069187108185302E+00 
-.7758670331396409E-h00 

7439849501 152674E-h00 
-.7114064976175457E4-00 
-.6782391686910609E-t-00 
-.6445701594098660E-f-00 
-.6104710013384929E-HOO 
-.5760010202980960E-HOO 
-.5412099413257457E+00 
-.5061398697742787E-h00 
-.4708268134473433E-h00 
-.4353018643598344E+00 

-  .3995921259242572E-H00 
-.3637214481257228E-F00 
-.3277110167114320E-1-00 
-.2915798305819667E+00 
-.2553450930388687E-KOO 
-.2190225363501577E-HOO 
-.1826266945721476E-1-00 
-.1461711362450572E4-00 
- .  1096686661347072E-h00 
-.7313150339365902E-01 
-.3657144220122915E-01 
0 


Weight 


0.4374483371752129E-02 
0.9842619236149078E-02 
0.1463518300250369E-01 
0.1862396111287527E-01 
0.2184988739217138E-01 
0.2442858670932862E-01 
0.2648864579258096E-01 
0.2814375940413615E-01 
0.2948528624795690E-01 
0.3058356160435090E-01 
0.3149181066633766E-01 
0.3225015506203403E-01 
0.3288893713079314E-01 
0.3343126421620424E-01 
0.3389488931551181E-01 
0.3429358206877410E-01 
0.3463812513892117E-01 
0.3493704033879884E-01 
0.3519712095895683E-01 
0.3542382499917732E-01 
0.3562156808557525E-01 
0.3579394352776868E-01 
0.3594388900778062E-01 
0.3607381381247460E-01 
0.3618569660385742E-01 
0.36281 16095737887E-01 
0.3636153393399723E-01 
0.3642789154364812E-01 
0.3648109393796617E-01 
0.3652181242257066E-01 
0.3655054982303338E-01 
0.3656765531685031E-01 
0.3657333451556860E-01 
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Fast  Mathematical  Algorithms  &  Hardware  Corporation 


Ronald  R.  Coifman 
Vladimir  Rokhlin 


June  21,  2002 

Defense  Technical  Information  Center/OCP 
8725  John  J.  Kingman  Road,  Suite  0944 
Fort  Belvoir,  VA  22060-6218 

Re:  Contract  #:  F49620-98-C-0051 

To  Whom  It  May  Concern; 

In  regard  to  the  above  referenced  contract,  due  to  typographical  errors,  would  you  please  correct 
the  SF298  forms  and  the  reports  to  read  Contract  #:  F49620-^-C-0051,  not  F49620-97-C-0051. 
I  apologize  for  any  inconvenience  this  may  have  caused. 

If  you  have  any  questions  please  feel  free  to  call  at  the  above  number  or  email  fran@fmah.com 


1020  Sherman  Avenue 
Hamden,  CT  06514 
USA  (203)  248-8212 
USA  (203)  287-8765FAX 


Sincerely  yours, 

Frah  Kearns,  Business  Manager 
F.M.A.&  H.  Corporation 


